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FOREWORD 


This  report  was  prepared  by  Lockheed  Missiles  and  Space  Company,  Inc., 

Palo  Alto  Research  Laboratories,  3251  Hanover  Street,  Palo  Alto,  California, 
in  partial  fulfillment  of  the  requirements  under  Contract  F33615-76-C-3105. 

The  effort  was  Initiated  under  Project  2307,  "Research  in  Flight  Vehicle 
Structures,"  Task  2307N102,  "Research  in  the  Behavior  of  Metallic  and 
Composite  Components  of  Air  Frame  Structures."  The  project  monitor  for 
the  contract  was  Dr.  Narendra  S.  Khot  of  the  Structures  and  Dynamics  Division 
(AFWAL/FIBRA). 

The  report  covers  work  under  the  contract  concerned  with  the  efficiency 
of  different  discretization  procedures.  The  technical  work  under  the  contract 
was  performed  during  the  period  June  1976  through  October  1980.  Review  reports 
were  submitted  in  October  1980  and  the  final  report  in  March  1981.  Results 
from  separate  but  related  efforts  are  included  for  completeness.  These  include 
the  description  of  an  element  configuration  developed  under  the  LMSC  independent 
research  program  and  a  summary  of  the  state-of-the-art  performed  under  contract 
with  AFOSR. 

The  other  reports  published  under  this  contract  are  "Imperfection  Sensi¬ 
tivity  of  Optimized  Structures,"  (AFWAL-TR-80-3128),  "Panel  Optimization  with 
Integrated  Software  (POIS),"  (AFWAL-TR-80-3073,  Vol  I  and  II),  "Design  of 
Composite  Material  Structures  for  Buckling,  An  Evaluation  of  State-of-the-Art," 
(AFWAL-TR-81-31 02) ,  "Supplementary  Studies  on  the  Sensitivity  of  Optimized 
Structures,"  (AFWAL-TR-8I-3013). 
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Section  I 


INTRODUCTION 

The  structural  response  to  a  given  environment  is  determined  by  the  differen¬ 
tial  equations  of  motion  of  deformable  bodies.  Numerical  solutions  of  such 
equations  are  discussed  in  the  following  for  the  case  in  which  all  derivatives 
with  respect  to  time  can  be  discarded  so  that  the  equations  degenerate  into  a 
set  of  static  equilibrium  equations.  Additional  questions  that  nay  arise  in  a 
dynamic  environment  are  not  considered.  We  pay  special  attention  to  the 
behavior  of  shell  structures.  For  such  structures,  the  mathematical  problem 
is  defined  by  a  set  of  non-linear  partial  differential  equations  of  elliptic 
type  or  by  an  equivalent  energy  principle.  Analytic  solutions  of  such  problems 
for  a  reasonably  large  class  of  structural  configurations  are  not  within  the 
realm  of  the  possible.  Consequently,  the  mathematical  problem  is  recast  into 
a  numerical  problem  for  solution  on  the  computer.  This  transformation  can  be 
applied  to  the  equilibrium  equations  or  directly  to  seme  energy  expression. 

The  output  from  the  computer  consists  of  a  sequence  of  numbers,  in  some  way 
representing  the  functions  satisfying  equilibrium  equations  and  boundary  condi¬ 
tions.  If  the  solution  is  represented  by  a  linear  superposition  of  a  set  of 
"basis  functions",  then  the  components  of  the  output  vector  consist  of  the 
coefficients  in  this  series.  This  is  the  case  if  we  use  the  Galerkin  or 
Rayleigh-Ritz  procedures,  if  we  use  the  finite  difference  or  finite  element 
procedures,  the  solution  function  is  represented  by  its  values  at  a  number  of 
discrete  locations  within  the  structure .  Because  these  discretized  methods 
are  readily  applied  in  a  computer  program  for  a  general  type  of  structure, 
they  have  been  gaining  popularity.  This  applies  in  particular  to  the  finite 
element  method.  The  discretized  methods  make  use  of  numerical  differentiation 
and  numerical  integration.  A  review  of  these  operations  here  is  intended  to 
serve  as  a  background  to  a  discussion  of  the  different  options  that  are 
availble  for  numerical  analysis  of  shell  structures. 


Section  II 


NUMERICAL  DIFFERENTIATION 


Numerical  differentiation  consists  of  the  replacanent  of  the  derivatives  of  a 
function  by  difference  quotients  (or  finite  difference  expressions) .  Such 
expressions  are  generally  based  on  local  polynomial  approximations.  A  truncated 
Taylor  series  can  be  used  for  this  purpose.  In  the  one-dirtensional  case 
(one  space  variable)  the  series  is  of  the  form 

2  n 

f(x)  =  f(0)  +  xf'(O)  +  |-f"(0)  +  .  .  .  .^f(n)  (0)  +  R  (1) 


where  the  prime  denotes  differentiation  with  respect  to  x.  The  remainder 
R  represents  the  sum  of  the  terms  that  were  excluded  when  the  series  was 
truncated  after  the  (iw-l)th  term.  It  can  be  shown  that  (see  Reference  1, 
for  example)  that 


xn+1 

R±Wi 


F 


where 


(2) 


F  =  Max  [f(n  +  1}  (£)];  0  <  £  <  x 

This  bound  on  the  truncation  error  is  useful  for  estimates  of  the  accuracy 
of  the  output  vector. 

Finite  difference  expressions  may  be  derived  at  a  number  of  control  points. 
First,  the  function  values  at  a  number  of  discrete  node  points  are  expressed 
in  terms  of  the  derivatives  at  the  control  points  by  use  of  a  suitably  truncated 
Taylor  series.  In  Equation  (1) ,  it  is  assumed  that  x  =  0  at  the  control  point. 
The  procedure  leads  to  a  set  of  linear  equations  and  through  solution  of  this 
system,  expressions  for  the  derivatives  at  the  control  point  are  obtained  in 
terms  of  the  function  values  at  the  selected  node  points.  The  number  of  node 
point  values  included  must  be  equal  to  the  number  of  terms  of  the  Taylor 
series.  In  the  general  case  the  highest  order  derivative  so  determined  will 


be  of  first  order  accuracy,  i.e.,  the  error  E  =  0(h) .  if  ex¬ 
pansion  is  complete  through  derivatives  of  nth  order,  the  kth  derivative  (k<n) 
will  be  of  order  (n  -  k  +  1) . 


Figure  1  illustrates  how  finite  difference  expressions  can  be  derived  in  a 
two-dimensional  space.  A  Taylor  series  approach  gives  us 
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Where  a  prime  indicates  differentiation  with  respect  to  the  coordinate, 
and  a  dot  differentiation  with  respect  to  the  coordinate. 


If  derivatives  up  to  the  second  order  are  to  be  determined,  it  is  sufficient 
to  specify  for  each  control  point  a  set  of  six  neighboring  node  points 
(Ni  to  Ni+5) .  By  applying  Equation  (3)  at  each  of  the  six  node  points,  we 
obtain  the  equation  system 
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D  Finite  Difference  Grid 


The  solution  of  the  equation  system  yields  a  set  of  finite  difference  expres¬ 
sions  for  the  derivatives  at  control  point  j.  The  error  bound  for  the  second 
order  derivatives  are  E  =  0(h) .  The  first  order  derivatives  are  of  second 
order  accuracy.  If  one  or  more  of  the  second  order  derivatives  is  left  out, 
the  remaining  derivatives  nay  be  expressed  in  terms  of  five  nodal  values. 

Then  the  first  order  derivatives  are  only  of  first  order  accuracy  and  the 
expressions  for  the  second  order  derivatives  may  be  meaningless,  E  =  O(h^) . 

Since  there  are  four  different  third  order  derivatives,  we  must  include  ten 

3 

nodal  function  values  in  order  to  raise  the  accuracy  by  one  order  [E  =  0(h  ) 

2 

for  the  first  and  E  =  0(h  )  for  second  order  derivaties] . 

Linear,  quadratic  or  cubic  expansions  are  obtained  if  all  derivatives  up  to 
first,  second  or  third  order  are  retained.  Bilinear,  biquadratic  or  bicubic 
are  expansions  including  all  derivatives  that  are  at  most  of  first,  second 
or  third  order  with  respect  to  any  one  of  the  two  space  variables.  It  may  be 
noticed  that  the  order  of  accuracy  of  bilinear  expansion  remains  the  same  as 
that  of  the  linear  expansion  and  a  similar  statement  holds  for  higher  order 
expansions. 

The  nodal  values  (degrees  of  freedom)  need  not  be  restricted  to  function 
values.  A  higher  order  derivative  at  a  control  point  can  be  expressed  in 
terms  of  lower  order  derivatives  at  the  nodal  points.  These  lower  order 
derivatives  are  then  included  as  support  for  the  local  approximation  of  the 
solution.  Such  procedures  are  discussed  in  Reference  2,  for  example. 

The  location  of  the  control  points  in  relation  to  the  node  points  can  be 
chosen  so  that  the  coefficient  becomes  zero  for  the  lowest  order  terms  among 
those  that  are  discarded.  In  that  case  the  accuracy  of  the  approximation 
is  raised  by  one  order.  For  demonstration  we  consider  the  configuration  in 
Figure  2. 

The  degrees  of  freedom  of  the  system  are  the  function  values  f1  and  f  2  and  the 
first  order  derivatives  and  at  the  node  points.  The  finite  difference 
equivalent  of  the  second  order  derivative  at  same  point  in  the  interval  is 


The  solution  of  this  equation  system  includes 


8 


between  two  nodes.  On  the  other  hand,  if  we  derive  a  formula  for  the  first 
derivative  at  a  node  point  fran  two  function  values  (forward  or  backward 
differences)  we  find  that  the  expression  is  of  first  order  accuracy,  i.e.. 


f.  if..  + 

h  2  1 
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By  use  of  three  node  points  we  can  define  the  first  order  derivative  with 
second  order  accuracy  at  any  point.  If  the  control  point  is  placed  at  the 
node  point  in  the  middle  the  derivative  is  defined  by  the  first  of 
Equations  (9) .  That  is,  the  error  is  four  times  larger  than  it  is  in  the 
derivative  at  the  stress  window  in  the  two  point  scheme.  On  the  other  hand 
the  forward  and  backward  difference  expressions  are  still  of  second  order 
accuracy 
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It  can  be  shown  that  first  order  derivatives  determined  from  a  uniform 
three  point  scheme  are  most  accurate  (of  third  order)  at  stress  windows, 
located  synmetrically  at  a  distance  of  h/2/3  from  the  nddooint  node. 


A  comparison  of  second  order  derivatives  frcm  three-  and  four-point  schemes 
gives  similar  results.  Ever  order  derivatives  in  a  uniform  grid  are  most 
accurately  computed  if  one  node  coincides  with  the  control  point  and  the 
remaining  nodes  are  placed  symmetrically  around  this  point.  The  location 
of  stress  windows  in  the  case  of  nonuniform  spacing  is  discussed  in 
Reference  3. 


Section  III 


NUMERICAL  INTEGRATION 

The  purpose  of  numerical  integration  is  to  compute  an  approximate  value  of 
b 

/  f(x)  dx  . 
a 

In  order  to  achieve  this  we  divide  the  range  of  the  integral  a,b  into  a 
number  of  small  subintervals.  If  the  function  values  are  determined  at  the 
midpoints  of  intervals  (see  Figure  5a) ,  the  integral  can  be  obtained  by 
use  of  the  so-called  rectangle  rule,  i.e. , 

b  N 

/  f (x)  dx  =  h  £  f.  (i: 

a  i=l  1 

where  N  is  the  number  of  intervals.  If  the  function  values  are  determined 
at  points  of  division  between  the  intervals  (including  the  end  points  a,b) , 
then  we  can  use  the  trapezoidal  rule  illustrated  in  Figure  5b. 


b 

/  f (x)  dx  =  h(V2  fQ  +  f.  +  f2  ...  f  x  +  1/2  f  ) 


It  is  shown  in  Reference  1  that  both  these  methods  are  of  second  order  accuracy. 
The  rectangle  rule  is  somewhat  mare  accurate  with  the  error  bound 


2 

E  <  1  ~  ~"24^~ '  I  f ' '  (5)|;  a  <  C  <  b. 
while  the  error  bound  for  the  trapezoidal  rule  is 
E  <  1  -  -  -^-2  |f"  (0|;  a  <  5  £  b. 


(15) 


(16) 


A  number  of  procedures  of  higher  accuracy  have  been  proposed  (Euler-McLaurin, 
Stirling) .  The  Newton-Cotes  series  of  integration  formulas  are  based  upon  the 
passing  of  a  polynomial  through  a  sequence  of  function  values  and  integration 
of  this  polynomial  over  the  subintervals.  Since  the  trapezoidal  rule  is 
based  on  a  linear  approximation  over  the  sub intervals ,  it  may  be  considered 
as  the  lowest  order  method  in  the  Newton-Cotes  series.  A  Newton-Cotes 
formula  of  third  order  accuracy  is  obtained  if  we  use  second  order  poly¬ 
nomials  for  interpolation  between  node  points.  With  uniform  spacing,  this 
member  of  the  Newton-Cotes  family  is  well  known  as  Simpson's  formula 


b  . 

/  f (x)  dx  =  ^  (fQ  +  4^  +  2f2  +  ...  4fn_x  +  fn)  (17) 

a 

We  notice  that  use  of  Simpson's  formula  requires  that  the  number  of  sub¬ 
intervals  is  even.  Higher  order  Newton-Cotes  formulas  are  increasingly 
restrictive  with  regard  to  the  permissible  number  of  subintervals. 


A  method  of  special  importance  in  finite  element  analysis  is  the  Gaussian 
Quadrature.  In  this  procedure  strategic  positions  are  established  for  the 
points  at  which  the  function  is  to  be  determined.  As  an  example  we  will 
show  here  hew  a  two-point  integration  scheme  is  derived.  If  the  function 
values  at  x  =  ±al  are  f^  and  f2  with  x  =  0  at  the  midpoint  of  the  interval, 
we  can  approximate  the  function  by 


f  = 


fl  +  f2 


x  +  F2(x)  +  R 


(18) 
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where  F2(x)  is  a  function  of  second  degree  in  x  and  vanishes  at 
x  =  ±  a£  and  R  is  a  power  series  beginning  with  the  third  order  term. 


Due  to  the  symmetry,  the  integral  of  odd  power  terms  vanishes  and  consequently 
we  can  write 


f  =  - 1--2  -  +  Cx  (o2£2  -  x2)  +  C2(a)x4  + 


(19) 


and 


/  f (x)  dx  =  4  (f.  +  f,)  +  C,£3  (a2  -  i)  +  C0  (a) £5/5  +  ...  (20) 

-i/2  *  1  £  i  12  * 

Choosing  to  determine  the  function  values  at  the  coordinates  corresponding 
to  a  =  ±  1/2/3  we  obtain 

1/2 

f  f(x)  dx  =  £(0.5f-,  +  0.5f_)  (1  +  E).  (21) 

-i/2 

4 

The  relative  error  E  =  C£  /(f^  +  f2)  where  C  is  a  constant.  The  method 
integrates  any  polynomial  up  to  the  third  order  exactly.  The  points  corre¬ 
sponding  to  a  =  ±1/2/5  are  referred  to  as  the  Gaussian  points.  Integration 
schemes  of  higher  order  are  obtained  by  inclusion  of  additional  Gaussian 
points.  Tables  are  available  that  give  the  location  of  the  Gaussian  points 
and  corresponding  weighting  factors  for  the  function  value.  For  exact 
integration  of  a  jth  order  polynomial  where  j  is  odd  we  need  (j  +  l)/2 
integration  points.  Sometimes  polynomials  are  integrated  with  fewer  points 
although  this  permits  the  existence  of  positive  definite  functions  whose 
integral  over  the  domain  vanishes .  The  procedure  is  then  referred  to  as 
reduced  integration. 

We  notice  that  the  position  of  the  "Gaussian  points"  coincides  with  the  stress 
windows  far  the  second  order  derivative  in  the  case  with  displacements  and 
rotations  defined  at  the  end  points  of  the  interval  (Figure  2) .  The  rectangle 


12 


rule,  identical  to  the  use  of  a  Gaussian  Quadrature  over  the  individual 
intervals  with  only  one  integration  point  in  each,  integrates  only  first  order 
functions  exactly. 


Section  IV 
ENERGY  METOCDS 

The  equations  of  equilibrium  of  a  deformable  body  can  be  derived  by  considera¬ 
tion  of  the  balance  between  elastic  forces,  body  forces  and  possible  surface 
tractions  in  an  infinitesimal  volume  element.  By  use  of  the  constitutive 
equations,  relating  strains  to  stresses,  and  the  kinematic  relations,  defining 
strains  in  terms  of  displacement  components,  the  equilibrium  equations  can  be 
expressed  in  terms  of  the  displacements. 

If  the  loading  is  conservative  and  if  only  elastic  deformation  is  considered, 
the  condition  stated  by  the  equilibrium  equations  is  equivalent  to  the  require¬ 
ment  that  the  sum  of  the  stored  strain  energy  and  the  potential  energy  of  the 
applied  force  system  equals  zero  during  a  snail  virtual  displacement.  Based 
on  this  theorem  of  virtual  work  we  can  derive  the  theorem  of  minimum  potential 
energy.  According  to  this  theorem: 

Among  admissible  displacements,  the  configuration 
corresponding  to  equilibrium  is  such  that  the  total 
potential  energy  is  stationary;  for  stable  equilibrium 
it  is  at  a  mininum. 

That  is 

5  (U  +  W)  =0  (22) 

where  U  represents  the  strain  energy  and  W  the  potential  energy  of  the 
external  forces,  i.e.  the  negative  of  the  work  dene  by  these  forces  during 
deformation.  Admissible  are  all  displacements  that  satisfy  continuity  and 
essential  boundary  conditions .  The  strain  energy  U  can  be  written 

U  =  1/2  /  e  [E]  c  dV  (23) 

V 
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where  c  represents  a  vector  of  strains  and  [E]  is  the  (3  x  3)  matrix  of 
coefficients  in  the  constitutive  relations.  The  potential  energy  of  the 
external  forces  is  represented  by 


W  =  -  /  U-T  ds  (24) 

ST 

where  T  represents  the  components  of  surface  traction  and  ST  the  part 
of  the  structure  on  which  such  tractions  (external  forces)  are  defined. 

The  theorem  of  minimum  complementary  energy  is  an  alternative  formulation 
also  derivable  from  the  principle  of  virtual  work.  This  theorem  states  that: 

Among  all  admissible  states  of  stress  the  configuration 
corresponding  to  equilibrium  is  such  that  the  complemen¬ 
tary  energy  is  stationary;  for  stable  equilibrium  it  is 
at  a  minimum. 

That  is 


6  (U*  +  W*)  =  0 


(25) 


where 


U*  =  1/2  J  a[E] -1a  dV 
V 


and 


W*  =  -  /  T'U  ds  (26) 
SU 

where  Sy  is  the  surface  upon  which  displacements  are  defined  and  T  the 
corresponding  tractions  (reactions) .  Admissible  are  stress  states  that 
satisfy  continuity  and  boundary  conditions  on  surface  tractions. 


Multifield  energy  principles  including  both  stresses  and  strains  as  free 
variables  can  be  derived  based  on  the  relation 

I  a  edV=U+U*  (27) 

V 

The  expression  for  the  energy  is  in  mathematical  terms  referred  to  as  a 
functional .  That  is,  it  defines  a  scalar  value  U  corresponding  to  a  given 
function.  In  analysis  of  plane  stress  in  plates  or  in  a  three  dimensional 
analysis,  the  functional  contains  derivatives  up  to  the  first  order  only. 

In  beam  analysis  or  in  analysis  of  plates  and  shells  certain  approximations 
result  in  the  introduction  of  second  order  derivatives  into  the  energy 
functional.  Displacements  at  any  point  in  the  structure  are  expressed  in 
terms  of  displacements  and  displacement  derivatives  at  a  reference  surface . 

As  a  consequence,  the  problem  is  reduaad  to  be  one-dimensional  (response 
quantities  are  functions  of  one  spatial  coordinate  only)  for  the  case  of 
beams  and  two-dimensional  for  plates  and  shells.  For  beams  then,  the  equa¬ 
tions  of  equilibrium  beacme  ordinary  differential  equations,  while  for  plates 
and  shells  they  retain  partial  differential  equations.  The  reduction  of  the 
number  of  spatial  coordinates  can  be  obtained  because  the  structure  is  thin, 
allowing  the  approximations: 

(1)  normals  to  the  reference  surface  remain  straight  and  normal 
during  deformation,  and 

(2)  the  transverse  normal  stress  is  negligibly  snail. 

In  a  second  order  theory,  approximately  accounting  for  transverse  shear 
deformation,  the  normals  remain  straight  but  are  allowed  to  rotate  in  rela¬ 
tion  to  the  reference  surface. 

The  assumptions  on  which  a  shell  theory  is  based  allows  the  integration  through 
the  thickness  so  that  the  strain  energy  is  represented  by  a  surface  integral. 
The  strain  energy  is  expressed  in  terms  of  a  six -component  vector  of  reference 
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surface  strains  and  changes  of  curvature  and  the  (6x6)  matrix  of  coefficients 
in  the  equations  relating  those  to  the  vector  of  force  and  moment  resultants. 


Solutions  to  the  equilibria  problem  are  obtained  by  application  of  one  of  the 
energy  theorems  after  numerical  analysis  procedures  has  been  used  to  express 
the  corresponding  functional  in  terms  of  a  finite  nunber  of  degrees  of 
freedom.  Another  possibility  is  to  seek  numerical  solutions  to  the  equivalent 
differential  equations  of  equilibrium. 


Section  V 


GLOBAL  FUNCTION  APPROACH 


In  this  section  we  briefly  discuss  the  methods  in  which  the  solution  function 
is  represented  by  a  linear  superposition  of  basis  functions.  These  are 
generally  allowed  to  assume  nonzero  values  over  the  entire  domain  (except  as 
restricted  by  boundary  conditions) .  We  refer  to  such  a  procedure  as  a 
global  function  approach  in  order  to  emphasize  its  distinction  from  the  finite 
element  method  to  be  discussed  later. 


The  Galerkin  Method  has  been  widely  used  to  yield  approximate  solutions  to 
differential  equations  in  structural  analysis  as  well  as  in  many  other 
applications.  The  method  is  applicable  to  partial  as  well  as  to  ordinary 
differential  equations.  For  simplicity,  we  write  the  equation  in  the  form 

L(u)  =  0  (28) 


where  L  is  a  differential  operator  and  u  represents  the  displacement  field. 
We  seek  a  solution  in  the  space  of  all  trial  functions  defined  by 

^=2.  Vn  (29) 

n=l 

where  the  basis  functions,  <t>n,  represent  kinematically  admissible  functions 
(i.e. ,  they  are  continuous  and  satisfy  given  boundary  conditions) .  The 
components  of  the  output  vector  are  the  an>  Applying  Galerkin' s  method, 
we  determine  the  an  through  solution  of  the  equation  system: 


/  L(u)  <f>  dV  =  0 
V  m 


<t>  dV  =  0  for  m  =  1,N. 
m 


(30) 
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If  u  =  u  represents  the  solution  of  this  system  with  N  basis  functions, 
convergence  to  the  correct  solution  is  implied  if  the  norm 


Hl^-u*!!  +  0  as  N  +  ®  (31) 

where  u*  is  the  solution  of  the  mathematical  problem  L(u)  =0.  The  norm  of  a 
function  defined  on  the  domain  V  can  be  chosen  as 

i |u| |  =  /  u2  dV  (32) 

V 

We  will  return  below  to  the  requirements  for  convergence  with  increasing 
value  of  N. 


The  energy  methods  can  also  be  used  directly  for  construction  of  a  solution 
to  the  euqilibrium  problem  in  terms  of  a  linear  combination  of  global 
functions.  One  such  procedure  is  the  Rayleigh-Ritz  Method  in  which  the 
trial  functions 


1?  = 


N 

£ 


a  * 
nyn 


(33) 


are  substituted  into  the  expression  for  the  total  potential  energy.  The 
basis  functions  <J>n  are  required  to  satisfy  essential  (displacement)  boundary 
conditions.  Hie  unknown  coefficients  a^  are  determined  through  minimization 
of  the  total  potential  energy.  Natural  boundary  conditions  are  automatically 
satisfied  through  the  minimization. 


A  proof  for  convergence  of  the  method  is  provided  in  Reference  4.  The  condi¬ 
tions  for  convergence  to  the  correct  solution  are 

(1)  For  functionals  containing  derivatives  up  to  the  nth  order  it  is 
required  that  derivatives  of  the  trial  functions  up  to  the  (n-l)th 
order  are  continuous.  This  is  a  sufficient,  but  as  we  will  see 
later,  not  always  necessary  continuity  condition. 
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(2)  The  set  of  trial  functions  must  be  complete.  That  is,  the 
set  of  trial  functions  must  contain  a  sequence  of  functions 
that  in  the  limit  approaches  any  admissible  function 
arbitrarily  close. 

(3)  Essential  (displacement)  boundary  conditions  are  satisfied. 

The  Galerkin  method  discussed  above  is  in  one  respect  an  extension  of  the 
Rayleigh-Ritz  method.  It  is  applied  also  to  differential  equations  that 
cannot  be  derived  through  a  variational  approach .  It  is  shown  in  Reference  4 
that  the  Galerkin  method  when  applied  to  variational  problems  with  quadratic 
functionals  is  identical  to  the  Rayleigh-Ritz  method.  This  establishes 
convergence  of  the  Galerkin  method  within  that  range.  For  other  cases, 
it  appears  that  the  assumption  of  convergence  is  based  on  conjecture. 

Before  the  introduction  of  the  digital  computer,  the  Rayleigh-Ritz  and  Galerkin 
procedures  were  frequently  applied  with  trial  functions  chosen  on  an  intuitive 
hasis.  Very  few  terms  were  included.  The  use  of  complete  sequences  and 
convergence  control  became  popular  with  the  arrival  of  high-speed  computers. 
With  the  demand  for  programs  of  general-purpose  type,  these  methods  were 
gradually  abandoned  in  favor  of  finite  difference  and  finite  element  methods. 
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SectiOi.  VI 


FINITE  DIFFERENCE  METHODS 

The  finite  difference  methods  are  based  on  replacement  of  derivatives  by  finite 
difference  expressions  as  discussed  above.  The  different  expressions  can  be 
introduced  into  the  equilibrium  equations  or  in  some  energy  expression.  In 
the  former  case  an  algebraic  equation  system  is  directly  obtained,  in  which 
the  node  point  displacements  and  possibly  rotations  are  the  unknowns  and  each 
equation  expresses  equilibrium  at  one  of  the  control  points.  The  number  of 
equations  must  be  equal  to  the  number  of  degrees  of  freedom  of  the  systan. 

Many  examples  of  application  of  this  procedure  are  given  in  the  literature. 

In  Reference  5,  it  is  applied  in  an  analysis  of  oolunn  buckling. 

The  problems  involved  in  practical  finite  difference  analysis  are  much  the 
same  if  finite  difference  expressions  are  used  in  combination  with  the  energy 
approach.  Since  this  approach  presently  is  of  more  direct  interest  to  us, 
the  discussion  of  the  direct  use  of  the  equilibrium  equations  will  be 
restricted  to  the  problem  of  convergence.  Over  the  years  many  efforts  have 
been  made  to  shew  that  the  method  converges  with  decreasing  node  point  spac¬ 
ing  to  the  solution  of  the  differential  equation,  i.e., 

Ml^-u*!!  -  Oash-vo  (34) 

where  U*1  is  the  finite  difference  solution  corresponding  to  a  node  point  spacing 
of  h. 

Rigorous  mathematical  proofs  have  been  presented  for  special  cases,  but 
due  to  the  diversity  of  differential  equation  forms,  boundary  conditions  and 
shapes  of  the  domain,  it  appears  to  be  difficult  to  establish  convergence 
in  the  general  case.  However,  mathematical  rigor  has  never  been  the  trade¬ 
mark  of  engineering  analysis.  Often  the  assumptions  and  simplifications  in 
modeling  the  structure  are  of  such  a  nature  as  to  make  the  quest  for  mathe¬ 
matical  rigor  rather  extravagant.  If  the  application  of  apparently  reasonable 
methods  had  been  deferred  in  anticipation  of  rigorous  proofs  of  uniqueness 


and  convergence,  progress  in  engineering  science  would  have  been  considerably 
retarded. 

The  finite  difference  approach  appears  reasonable  from  the  point  of  view  of  an 
engineering  analyst.  One  difficulty  related  to  the  convergence  proof  is  that 
the  error  term  in  the  Taylor  series  expansion  contains  a  derivative  [fn+^‘]  U) 
in  Equation  2  of  the  solution  function  itself.  Witn  seme  feeling  for  the 
physical  behavior  of  the  system,  the  analyst  nay  assume  that  the  solution 
vector  varies  continuously  with  the  input  data  ( |fn+1,  (0  |<  C  j  |T|  |)  so  that 
the  truncation  error  can  be  expressed  in  terms  of  the  loading  function  rather 
than  in  terms  of  the  solution.  Whenever  the  assumption  that  the  solution 
varies  continuously  with  the  input  data  is  violated,  finite  difference  formula¬ 
tions  may  lead  to  spurious  results.  Consider,  for  example,  the  case  of  a 
beam  in  bending.  We  substitute  the  second  order  derivative  by  a  three-point 
central  finite  difference  expression  except  at  one  internal  point  where  we 
use  either  a  backward  or  a  forward  difference  expression.  The  beam  so  defined 
has  a  link  at  the  point  in  question. 

It  appears  reasonable  to  accept  the  proposition  in  Reference  6  that  the  finite 
difference  approach  in  solution  of  differential  equations  converges  toward 
the  correct  solution  if: 

(1)  the  local  truncation  error  vanishes  with  the  grid  size 

(2)  for  small  values  of  h  the  solution  varies  continuously 
with  the  input  data  (loads) . 

The  second  requirement  will  exclude  a  spurious  solution  such  as  the  one  for 
the  beam  discussed  above. 

The  use  of  finite  differences  in  a  variational  approach  is  discussed  in 
Reference  7  on  page  182.  The  finite  difference  expressions  are  introduced 
directly  into  the  energy  expression  and  the  potential  energy  is  minimized 
with  respect  to  the  values  of  nodal  displacement  components .  The  convergence 
of  the  procedure  is  not  discussed  in  Reference  7.  We  nay  appeal  to  equiva¬ 
lence  with  finite  difference  solution  of  the  Euler  equations  discussed  above. 
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A  rigorous  and  general  proof  of  convergence  does  not  appear  bo  be  available. 

It  will  suffice  for  us  that  whenever  the  procedure  has  been  applied  to  a 
case  with  a  known  solution,  it  has  been  found  to  converge  toward  that  solution. 


The  advantages  of  introducing  the  difference  quotients  into  the  potential 
energy  expression  rather  than  into  the  equilibrium  equations  cure  that  the 
coefficient  matrix  becomes  synmetric  and  that  the  natural  boundary  conditions 
are  automatically  satisfied.  Ftor  a  simple  demonstration  of  the  method,  we 
will  consider  the  buckling  of  a  oolunn  with  uniform  cross-section  (compare 
Reference  5,  page  283) .  The  buckling  load  is  defined  as  the  load  level  at 
which  the  second  variation  of  the  total  potential  energy  vanishes.  Hence, 


262V  -  S2  j/  (EI(w/xx)2  -  P(w,x)2]  dx  »  0  j 


(35) 


In  order  to  take  advantage  of  the  stress  windows  as  discussed  above  [see 
Equations  (9)  and  (10)  ] ,  we  will  determine  w,x  at  midnodes  and  w#xx  at  node 
points  as  shewn  in  Figure  6. 


X  X  X  X  X 


X  NODE  POINTS 
CONTROL  PTS  FOR 

O  CONTROL  PTS  FOR  iJyx 


Figure  6.  Finite  Difference  Scheme  far  Colunn  Buckling  Analysis 


2  2 
Then  (wlj{)  can  be  integrated  by  use  of  the  rectangle  rule,  and  (w,^)  by 

use  of  the  trapezoidal  rule.  We  have 


N-l 


J(w  )2dx  =  h£  F(wi+1  -  wi} 
i=l 


(36) 


and 


N-l 


/(w,  )  2dx  =  ^  (w2  “  2wx  +  wQ)  2  +  h  £  |  (wi+1  -  2wi  +  w._1) 

1=2 


+  2h  (WN-1  “  ^N  +  WN+1) 


(37) 


Fictitious  points  (corresponding  bo  wQ  and  wN+1)  have  been  introduced  so  that 
the  second  order  derivatives  can  be  determined  at  the  end  points.  We  could 
instead  introduce  forward  and  backward  derivatives  at  these  points  or  add 
the  rotation  at  the  end  points  as  a  freedom.  Whenever  energy  methods  are  used, 
the  control  points  at  which  derivatives  are  required  coincide  with  the  integra¬ 
tion  points  at  which  the  energy  density  is  defined. 


After  the  number  of  uniformly  spaced  node  points  have  been  chosen,  we  can 
form  a  homogeneous  equation  system  in  which  P  appears  as  the  eigenvalue 
parameter.  Solutions,  given  in  Reference  5,  are  shewn  in  Table  1. 


TABLE  1 

COLUMN  BUCKLING  LOADS 


Number  of  Nodes 
on  Half  Column  (i) 


PCR 

Critical  Load  /  [EI(VL)  '1 


3 

4 

5 
7 
9 

11 


0.9495 

0.9774 

0.9372 

0.9943 

0.9968 

0.9979 

1.0000 
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Since  the  finite  difference  expressions  as  well  as  the  numerical  integrations 
are  of  second  order  accuracy,  we  expect  the  error  in  the  solution  to  be 
proportional  to  the  square  of  the  spacing  between  node  points  if  the  grid  is 
fine  enough  so  that  fourth  order  terms  in  the  error  are  insignificant. 

(We  have  emitted  the  step  shewing  that  natural  boundary  conditions  are  satisfied 
to  the  same  degree  of  accuracy.)  As  the  grid  spacing  equals  1/2 (i— 1)  ,  we  obtain 
frem  the  first  two  results  in  the  table  above 


PCR  -  °*9495 
PCR  “  °-9774 


(38) 


Dividing  the  first  of  Eqations  (38)  by  the  second,  we  obtain  an  equation  from 

which  we  determine  P  =  0.9997.  in  view  of  the  fact  that  the  values  of  P_ 

ck  CR 

for  3  or  4  points  were  rounded  to  four  figures,  this  is  as  close  to  the  exact 
solution  as  we  possibly  can  expect  and  thus  verifies  the  assumption  of  second 
order  convergence.  First  order  derivatives  frem  two  nodal  points  and  second 
order  derivatives  from  three  are  generally  of  first  order  accuracy.  A  second 
order  accuracy  is  obtained  by  a  favorable  choice  of  nodes  and  integration 
points.  This  type  of  phenomenon  is  referred  to  as  superoonvergence  in 
Reference  8.  By  use  of  two  solutions  corresponding  to  very  coarse  nodal 
spacing,  and  the  assumption  of  second  order  accuracy,  we  are  in  this  case 
able  to  predict  a  very  accurate  result  through  extrapolation.  This  method 
is  referred  to  as  Richardson's  extrapolation. 


In  the  two-dimensional  case,  it  is  more  difficult  to  utilize  stress  windows 
in  order  to  obtain  superconvergence.  We  will  briefly  discuss  the  problems 
involved  in  the  definition  of  efficient  finite  difference  schemes  for 
rectangular  nets  with  uniform  spacing.  (In  the  case  of  nonuniform  spacing, 
it  is  possible  to  make  use  of  the  same  expressions  by  mapping  the  shell  upon  a 
dorain  with  a  suitable  definition  of  a  distance.)  Examples  of  two-dimensional 
finite  difference  schemes  ctre  shewn  in  Figure  7.  The  scheme  used  in  Figure  7a 
defines  the  inplane  displacements  at  half-stations.  A  rectangular  integration 
scheme  is  used,  and  the  strain  energy  is  defined  at  the  w-node  at  the  center 


25 


SCHEME  1  (Half -Station) 


x  w  Defined 
o  u,v  Defined 


Integration  Area 


SCHEME  2  (Half -Station  by  NooR.) 
-x - e - x — »y,v 

x  w  Defined 
o  u  Defined 
□  v  Defined 


Integration  Area  For  Shear  u  ,  v  , w 


,*y 


Integration  Area  For  Direct  (u  .,vulw„,  etc.) 

,x  ijr  |*A 


SCHEME  3  (Whole-Station) 


•XvXvX;:*: 

x  u,v,w  Defined 

Integration  Area  For  Membrane  Strain 
Energy  u  x,  vy,uty,wx  ,  etc. 

'Integration  Area  For  Bending  Strain 

Energy 


Figure  7.  TVro-D  Finite  Difference  Schemes 
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of  the  figure.  This  point  then  represents  a  stress  window  for  bending  as  well 
as  membrane  strains.  There  is  one  obvious  weakness  in  this  s<  v'  me.  The 
membrane  strain  e^,  for  example,  at  the  integration  point  must  be  obtained 
as  an  average  of  its  values  at  points  A  and  B.  .While  this  operation  is 
still  of  second  order  accuracy  the  coefficient  of  the  error  term  is  increased. 

This  problem  is  eliminated  in  the  modified  half-station  scheme  introduced  in 
Reference  9  by  Noor  (Figure  7b) .  However,  both  the  half-station  schemes  shown 
here  are  somewhat  inefficient  for  nonlinear  or  stability  analysis.  It  is 
sufficient  for  illustration  of  this  problem  to  consider  a  straight  beam  element 


ex  -  u'x  +  V2(u?x  ♦  «?x> 


For  definition  of  the  membrane  strain  ex,  we  need  to  express  both  the  spatial 
derivatives  u,  and  w,  at  the  same  integration  points.  If  w  is  defined  at 
the  node  points  and  u  at  half -stations,  it  is  not  possible  to  use  the  most 
favorable  expression  for  u,  as  well  as  for  w,  . 

By  introduction  of  the  whole-station,  illustrated  in  Figure  7c,  both  these 
problems  are  eliminated.  With  the  membrane  and  bending  energies  integrated 
over  different  sets  of  integration  points,  it  is  possible  to  make  use  of  the 
stress  window  for  u.  as  well  as  for  w,  .  This  results  in  much  better 

X  X 

convergence  for  the  buckling  load  as  illustrated  in  Table  2. 

TABLE  2 

BUCKLING  LOADS  VERSUS  GRID  SIZE  FOR 
AXIALLY  LOADED  CYLINDRICAL  PANELS 


Grid 

Axial  x  Circumf. 


10  x  5 
12  x  6 
16  x  8 
20  x  10 
32  x  16 


Half -Station 
Figure  7a 


Whole  Station 
Figure  7c 


We  conclude  from  the  results  in  Table  2  that  for  similar  accuracy  we  can  use 
almost  twice  as  large  grid  spacing  with  the  whole-station  scheme.  Also,  it 
appears  that  the  error  in  both  cases  varies  quadrat ically  with  the  grid  spacing. 
However,  the  eigenvalues  were  not  computed  with  sufficient  accuracy  to  give 
reliable  information  about  the  rate  of  convergence. 

Unfortunately,  the  whole-station  schane  becomes  less  efficient  if  the  constitu¬ 
tive  relations  contain  terms  that  couple  force  and  moment  resultants  (such  as 
eccentrically  stiffened  shells) .  This  problem  occurs  because  the  strain  energy 
then  includes  products  of  first  order  derivatives  of  in-plane  and  normal 
displacements.  It  appears  possible  to  circumvent  this  problem,  but  since  the 
interest  largely  has  shifted  to  finite  element  formulations,  this  search  has 
not  been  vigorously  pursued. 

Another  problem  is  related  to  the  integration  of  the  membrane  strain  energy 
by  use  of  the  rectangle  rule,  or  as  observed  above,  by  Gaussian  integration 
with  only  one  point.  A  polynomial  that  matches  the  function  values  at  all 
four  corners  must  include  at  least  one  second  order  term.  Therefore,  this  is 
a  case  of  reduced  integration  allowing  a  nonzero  deformation  pattern  with 
vanishing  strain  energy.  This  pattern  is  defined  by  u  =  c^xy  and  v  =  c2xy, 
where  x  and  y  are  spaoe  coordinates  with  x=y=0  at  the  integration  point. 
Such  a  deformation  is  referred  to  as  a  mechanism,  and  if  allowed,  sometimes 
leads  to  spurious  solutions.  The  deformation  pattern  is  usually  prevented 
by  displacement  constraints  on  shell  boundaries  and  therefore  it  is  somewhat 
questionable  whether  it  is  worthwhile  to  remedy  the  situation  by  use  of 
integration  with  four  Gaussian  points. 


Section  vn 


FINITE  ELEMENT  ANALYSIS 


In  the  previous  paragraph  it  is  demonstrated  that  for  maximum  efficiency  in  a 
finite  difference  discretization  it  is  important  to  choose  the  position  of 
each  control  point  (or  integration  point)  judiciously.  It  is  important  also 
that  the  functions  and  their  derivatives  are  expressed  in  terms  of  the  function 
values  at  a  suitable  set  of  neighboring  nodes.  The  major  flaw  in  the  whole- 
station  scheme  of  the  previous  section  was  that  both  factors  in  terms 
representing  coupling  between  bending  and  membrane  action  cannot  be  represented 
by  the  most  efficient  expression  at  the  same  integration  points. 

Continuing  tht  search  for  a  better  formulation  we  may  consider  the  scheme 
illustrated  in  the  one  dimensional  case  in  Figure  8.  One  integration  interval 
is  shown  in  the  figure.  According  to  Equation  (8)  the  second  order  derivative 
of  the  lateral  displacement  is  of  third  order  accuracy  at  the  integration 
points.  The  first  order  derivative  of  the  lateral  displacement  based  on  the 
same  four  freedoms  is  at  least  of  third  order  accuracy  anywhere  in  the 
interval.  Based  on  three  function  values,  at  the  midpoint  and  at  the  end 
points  of  the  interval,  first  order  derivatives  (of  inplane  displacements) 
are  of  third  order  accuracy  at  the  Gaussian  points  for  two-point  integration 
[compare  the  discussion  following  Equation  (11) ] .  As  a  consequence  all 
quantities  included  in  the  strain  energy  are  of  third  order  accuracy  at  the 
two  Gaussian  points. 
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Integration  Area 


1/2 - 

M/(2/3  )- 


-{a— 6 


»l/(2/5)-H 


□  Integration  points 

X  Node  points  for  lateral  displacements 
and  their  first  order  derivatives 

O  Node  points  for  inplane  displacements 


Figure  8.  A  Finite  Difference  Scheme  with  Third  Order  Accuracy 

Characteristic  for  the  scheme  in  Figure  8  is  that  the  derivatives  at  all 
integration  points  within  any  integration  interval  are  based  on  the  same 
system  of  neighboring  nodal  values,  and  that  no  nodal  values  outside  of  the 
integration  area  are  engaged  in  the  formulation.  Due  to  these  properties  the 
formulation  adheres  in  the  strictest  sense  to  the  rules  that  define  a 
finite  element  procedure.  However,  these  rules  are  somewhat  artificial  and 
there  appears  to  be  no  meaningful  distinction  between  the  finite  element  and 
the  energy  based  finite  difference  analyses.  In  Reference  10  Felippa  refers 
to  this  finite  difference  procedure  as  finite  elements  with  extended  support 
because  the  energy  density  at  the  integration  points  is  expressed  in  terms  of 
nodal  displacement  freedcmes  outside  of  the  closed  domain  of  the  element. 

We  may  refer  to  configurations  that  satisfy  the  stricter  rule  as  self-contained 
elements . 

It  may  be  noticed  at  this  point  that  the  finite  difference  formulation  for  the 
membrane  energy  in  the  two-dimensional  scheme  is  a  self-contained  finite 
element.  In  the  bending  part  on  the  other  hand  the  energy  expression  is  based 
on  freedoms  at  nodes  outside  the  element  boundaries. 
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It  is  clear  then,  that  the  finite  element  method  could  have  been  derived 
through  a  specialization  or  refinement  of  the  finite  difference  energy  method. 

As  will  be  seen  in  the  fol.To.wing  it  could  also  have  been  the  result  of  a 
specialization  of  the  Rayleigh-Ritz  method,  in  either  case,  it  would  have 
been  directly  applicable  for  solution  of  any  set  of  ordinary  or  partial 
differential  equations  derivable  from  a  potential.  However,  in  contrast  to 
the  Rayleigh-Ritz  procedure  and  the  finite  difference  method  the  finite  element 
method  was  originally  derived  by  means  of  physical  and  largely  heuristic 
considerations  in  the  field  of  structural  mechanics  rather  than  in  the  field 
of  applied  m  there  tics.  As  an  afterthought  the  method  was  given  a  mathem¬ 
atical  interpretation  (Reference  11) .  This  led  to  considerable  refinement 
of  the  method  and  made  possible  the  extension  to  problems  outside  of  the 
field  of  structural  mechanics. 

The  mathematical  interpretation  of  the  finite  element  method  is  based  on  the 
Ritz  principle.  While  the  heuristic  mechanical  approach  may  have  great  appeal 
to  many  engineers,  the  best  understanding  of  the  method,  its  scope  and  its 
convergence  properties  is  obtained  if  it  is  presented  as  a  special  form  of  the 
Rayleigh-Ritz  method,  mich  in  the  same  way  as  in  the  book  by  Strang  and  Fix 
(Reference  6)  .  From  this  point  of  view,  the  finite  element  method  entails  the 
definition  of  a  set  of  trial  functions,  complete  in  the  sense  that  it  contains 
at  least  one  function  arbitrarily  close  to  any  admissible  solution  (displacement 
field) .  The  solution  of  the  problem  is  represented  by  that  member  of  the  set 
of  trial  functions  which  renders  the  functional  (energy  expression)  stationary. 
Typical  for  the  finite  element  method  is  that  the  trial  functions  are  obtained 
as  a  linear  combination  of  locally  defined  basis  functions .  That  is ,  each 
basis  function  is  zero  over  the  major  part  of  the  domain.  The  advantage  with 
this  arrangement  is  that  most  of  the  basis  functions  are  uncoupled,  and  thus, 
the  coefficient  matrix  of  the  final  equation  system  has  a  relatively  narrow 
bandwidth.  Also,  the  system  is  generally  well  conditioned. 
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Section  VIII 


ONE-DIMENSIONAL  PROBLEMS 


For  the  purpose  of  defining  local  basis  functions,  the  domain  is  subdivided 
into  a  number  of  intervals  or  subdomains  called  finite  elements.  Once  the 
domain  has  beer  so  divided,  the  unknown  function  within  each  elonent  can  be 
expressed,  for  ample,  in  terms  of  local  power  series  with  unknown  coefficients. 
Each  local  power  series  is  constrained  by  the  condition  that  at  the  endpoint 
of  the  interval  its  sum  equals  the  discrete  value  of  the  solution  functions 
f  ^  at  corresponding  node  point.  The  nodal  values  of  the  function  and  the 
free  coefficients  in  the  power  series  are  the  degrees  of  freedom  of  the  systan 
with  respect  bo  which  the  functional  is  minimized,  if  additional  constraints 
on  these  freedom  are  introduced  as  dictated  by  boundary  conditions,  the  trial 
functions  so  obtained  form  a  complete  set  of  admissible  functions.  With  a 
fixed  subdivision,  a  mesh,  the  finite  element  method  properly  applied  would 
converge  with  increasing  order  of  the  local  approximations.  In  such  a  case 
we  have  p-oonvergence .  If  the  mesh  is  coarse,  a  large  number  of  terms  may  be 
needed  in  the  power  series,  then  the  band  width  of  the  equation  system  becon.es 
large  and  also  the  system  may  become  ill-conditioned. 

A  more  cannon  practice  is  to  use  a  fixed  form  of  the  local  shape  functions . 

In  that  case  h-conyergence  is  obtained  with  gradually  refined  mesh  spacing. 

Or  rather,  the  mesh  is  made  fine  enough  so  the  analyst  feels  reasonably 
certain  of  obtaining  a  solution  of  satisfactory  accuracy. 

In  the  simple  one-dimensional  case  with  only  first  order  derivatives  the 
basis  functions  may  be  linear  as  shown  in  Figure  9.  The  value  of  $n  is  equal 
to  unity  at  node  n  and  zero  at  all  other  nodes. 


It  is  clear  then,  that  the  finite  element  method  oould  have  been  derived 
through  a  specialization  or  refinement  of  the  finite  difference  energy  method. 

As  will  be  seen  in  the  following  it  oould  also  have  been  the  result  of  a 
specialization  of  the  Rayleigh-Ritz  method.  In  either  case,  it  would  have 
been  directly  applicable  for  solution  of  any  set  of  ordinary  or  partial 
differential  equations  derivable  from  a  potential.  However,  in  contrast  to 
the  Rayleigh-Ritz  procedure  and  the  finite  difference  method  the  f  jiit e  element 
method  was  originally  derived  by  means  of  physical  and  largely  heuristic 
considerations  in  the  field  of  structural  mechanics  rather  than  in  the  field 
of  applied  mathenatics.  As  an  afterthought  the  method  was  given  a  mathen- 
atical  interpretation  (Reference  11) .  This  led  to  considerable  refinement 
of  the  method  and  made  possible  the  extension  to  problems  outside  of  the 
field  of  structural  mechanics. 

The  mathematical  interpretation  of  the  finite  element  method  is  based  on  the 
Ritz  principle.  While  the  heuristic  mechanical  approach  may  have  great  appeal 
to  many  engineers,  the  best  understanding  of  the  method,  its  scope  and  its 
convergence  properties  is  obtained  if  it  is  presented  as  a  special  form  of  the 
Rayleigh-Ritz  method,  much  in  the  same  way  as  in  the  book  by  Strang  and  Fix 
(Reference  6)  .  From  this  point  of  view,  the  finite  element  method  entails  the 
definition  of  a  set  of  trial  functions,  complete  in  the  sense  that  it  contains 
at  least  one  function  arbitrarily  close  to  any  admissible  solution  (displacement 
field) .  The  solution  of  the  problem  is  represented  by  that  msnber  of  the  set 
of  trial  functions  which  renders  the  functional  (energy  expression)  stationary . 
Typical  for  the  finite  element  method  is  that  the  trial  functions  are  obtained 
as  a  linear  ccmbiration  of  locally  defined  basis  functions.  That  is,  each 
basis  function  is  zero  over  the  major  part  of  the  domain.  The  advantage  with 
this  arrangement  is  that  most  of  the  basis  functions  are  uncoupled,  and  thus, 
the  coefficient  matrix  of  the  final  equation  system  has  a  relatively  narrow 
bandwidth.  Also,  the  systan  is  generally  well  conditioned. 
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Figure  9.  Linear  Basis  Function 


The  trial  functions  (with  f  =  =  0)  have  the  form 


„  N  -  1 

f  =  V*  a  <fc  where 
_  n  n 
n  =  2 


(40) 


(x  -  xn  -  l)/(xn  ■xn-l,ta\.li"S<n 
>n  -{  (xn  *  1  -  x)/(xn  *1  -  V  £ar  *n  <  x  i’S,  +  1 


0  elsewhere 


This  set  consists  of  all  piecewise  linear  functions.  It  is  complete  since, 
with  decreasing  mesh  spacing  (increasing  N)  any  continuous  function  can  be 
arbitrarily  closely  approximated.  The  trial  functions  (but  not  its  first 
derivatives)  are  continuous  over  the  domain  as  required  by  the  Ritz  theory. 

The  local  power  series  are  of  the  form 

f. -  f. 

f  =  fi  +  “  <*  -  V  (41) 

1  xi+l  -  x.  1 


That  is  they  are  of  first  order  in  x. 


The  trial  functions  defined  in  Eq.  (40)  may  be  used  to  determine  arbitrarily 
closely  the  function  f  (x)  that  minimizes  the  functional 

U  =  /  F(x)f?  +  G(x)f  dx  (42) 

L  x 

Where  L  is  the  total  length  of  the  domain  on  which  f  is  defined  and  F(x)  , 

G(x)  are  known  functions  x.  After  substitution  of  Equation  (40)  into 
Equation  (42)  the  integral  is  evaluated  over  each  interval  and  surrmed  over 
the  domain.  Due  to  the  local  characters  of  the  basis  functions,  the  products 
W  31:6  zero  everywhere  if  |n-m|<  1.  Consequently  the  terms  0^^  in  the 
functional  corresponding  to  F(x)f ,  will  vanish  if |n  -  m| >  1,  i.e.,  the 
equation  system  obtained  through  minimization  (3U/3an  =  0,  n  =  2,  N-l)  will 
be  narrowly  banded  (tridiagonal) . 

In  analysis  of  beams,  plates  and  shells  the  strain  energy  functional  includes 
second  order  derivatives  of  the  lateral  displacements.  For  such  problems  the 
Ritz  theory  requires  that  the  trial  functions  belong  to  the  function  space  C  . 

That  is,  it  is  required  that  the  trial  functions,  as  well  as  their  first  order 
derivatives,  are  continuous  over  the  domain.  In  the  one  dimensional  case 
(beam  analysis)  this  is  readily  achieved  if  the  function  value  as  well  as  its 
first  order  derivative  (rotation  in  beam  analysis)  are  considered  as  nodal 
freedoms.  A  cubic  representation  of  the  lateral  displacement  field  can  be 
obtained  by  use  of  the  Taylor  series  approach  [compare  Equations  (6)  and  (7)  ] . 

Continuity  requiranents  are  satisfied  since  the  values  of  the  functions  and 
their  first  order  derivatives  at  any  node  are  common  to  the  two  elements 
connected  at  the  node. 
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An  equivalent  representation  can  be  obtained  by  use  of  a  set  of  four  shape 
functions  defined  (Figure  10)  such  that  in  each  function  one  nodal  freedom 
equals  unity  and  the  other  three  are  zero.  The  basis  function  corresponding 
to  the  nth  node  is 

£n  ’  £  <Vn  +  bn*n>  <« 

where  a  is  the  function  value  and  b  its  first  order  derivative  at  node  n. 
n  n 

Here  an  and  bn  represent  the  degrees  of  freedom  of  the  system.  The  bending 
moments  and  the  rotations  at  each  of  the  integration  points  can  be  expressed 
in  terms  of  these  freedoms.  The  total  bending  strain  energy  is  then  readily 
expressed  in  the  form 

V  =  {g}T  [KK  {g}  (44 

e  c 

where 

«J>T-  “Vi  '  (ar>>2  ■  (bn>l  •  (bn>2’  (45 

and  K  is  referred  to  as  the  element  stiffness  matrix.  Summation  over  all 
elements  gives  sin  assembled  or  global  stiffness  matrix  (see  Reference  12, 
for  example) . 

Column  buckling  analysis  includes  such  terms  as  w,^  [see  Equation  (35)]. 

With  a  cubic  representation  of  w  this  term  is  a  fourth  order  polynomial, 
accurately  integrated  by  use  of  a  Gaussian  scheme  with  three  integration 
points. 

2 

The  critical  load  for  a  cantilevered  column  is  1/4  El  U/L)  .  A  finite 
element  buckling  analysis  of  a  column  with  (EI)/(4L^)  =  1.0  based  on  cubic 
representation  of  the  lateral  displacement  and  three  Gaussian  points  gives  the 
results  shown  in  Table  3. 


- - h 


Figure  10. 
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TABLE  3 

COLUMN  BUCKLING  RESULTS 


Number  of 
Elements,  i 

PCR 

ERROR 

E 

Ei4 

1 

9.943847 

.07424 

.0742 

2 

9.874659 

.005055 

.0809 

3 

9.870620 

.001016 

.0823 

4 

9.869928 

.000323 

.0827 

9 

9.869617 

.000013 

.0853 

2 

With  increasing  nunber  of  elements  the  accurate  solution,  P,_  =  ir  =  9.869604 

UK 

is  rapidly  approached.  It  appears  that  the  error  times  the  fourth  power  of 
the  number  of  elements  is  almost  a  constant.  This  would  indicate  that  the 
error  is  of  the  fourth  order  rather  than  third  as  expected,  may  be  somewhat 
fortuitously.  With  very  coarse  spacing  the  lowest  order  term  is  not  the 
dominating  error  and  with  very  fine  spacing  the  accuracy  may  be  governed  by 
round  off  errors.  As  a  compromise,  we  base  an  extrapolation  on  the  values 
corresponding  to  3  and  4  elements.  The  assumption  of  fourth  order  accuracy 
then  leads  bo  a  value  of  9.869612  for  the  critical  load,  i.e.,  the  relative 
error  is  less  than  10-^. 


With  only  one  element,  corresponding  to  two  elements  per  half-wave  in  the 
buckling  pattern,  the  error  is  still  less  than  one  percent.  In  engineering 
analysis  such  accuracy  is  generally  quite  acceptable.  However,  with  a  solution 
available  only  for  one  grid  size  the  analyst  does  not  knew  whether  it  is 
within  the  range  of  acceptable  accuracy.  Unless  a  sequence  of  gradually 
refined  solutions  is  available  he  is  not  on  firm  ground. 

For  analysis  of  rings  or  arches  it  seems  natural  to  develop  curved  beam 
elements.  In  that  case,  the  special  problem  of  strain  energy  due  to  a  rigid 
body  displacement,  self -straining,  must  be  considered.  For  illustration,  an 
element  of  a  circular  arch  is  shewn  in  Figure  11. 


If  the  circular  arch  is  rigidly  displaced  a  distance  6  in  the  x-direction, 
the  displacement  components  are  defined  by 


w  =  <S  sin  0 
v  =  6  cos  © 


(46) 


It  is  not  possible  with  a  truncated  power  series  to  represent  this  displace¬ 
ment  pattern  exactly.  As  a  consequence  a  rigid  body  displacement  introduces 
same  strain  energy  in  the  element.  This  situation  occurs  whenever  the 
geometry  of  the  element  cannot  be  exactly  represented  by  use  of  the  functions 
(usually  polynomials)  representing  the  displacement  components . 


Figure  11.  Rigid  Displacement  of  Curved  Element 

Elements  with  the  strain  energy  for  rigid  body  displacements  proportional  to 
seme  power  of  h  may  still  give  good  convergence  in  most  applications.  However 
there  are  cases  in  which  the  rigid  body  displacement  of  an  element  is  very  large 
in  comparison  to  the  displacements  corresponding  to  element  distortion.  In 
such  cases,  h-convergence  may  be  very  slow  for  elements  in  which  the  rigid 
body  energy  is  not  exactly  zero  but  proportional  to  sane  power  of  the  nodal 
spacing. 


Flat  elements  are  often  used  to  represent  curved  beams  or  surfaces.  The  reason 
for  this  approach  is  that  it  simplifies  the  formulation  and  also  that  it 
eliminates  problems  with  strain  energy  under  rigid  body  displacement.  Table  4 
shews  seme  results  from  a  buckling  analysis  with  cubic  straight  beam  elements 
of  a  ring  under  constant  direction  pressure.  The  critical  pressure  is 

4EI 

PCH  =  ^  (47) 
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The  results  in  Table  4  correspond  to  a  ring  with  EI/R  =0.25  lbs/in,  i.e. , 
the  exact  solution  is  p  =  1.0. 
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TABLE  4 

RING  BUCKLING  RESULTS 


Number  of 
elements ,  i 
(Over  90°  arc) 

PCR 

Error 

E 

Ei2 

4 

1.006970 

.006970 

.1015 

6 

1.002965 

.002968 

.1067 

8 

1.001641 

.001641 

.1050 

10 

1.001042 

.001042 

.1042 

We  notice  that  the  method  now  is  of  second  accuracy  (Ei  is  constant)  and  that 
to  achieve  somewhat  less  than  one  percent  error  it  is  necessary  to  use  four 
elements  per  half-wave  in  the  buckle  pattern.  Use  of  Richardson's  extrapola¬ 
tion  based  on  second  order  accuracy  with  8  and  10  elements  leads  p_,  =  1.00002 

vJR 

lbs/ in. 

As  long  as  each  element  is  self-contained,  the  order  of  the  accuracy  is  not 
changed  if  the  mesh  has  a  variable  spacing.  For  buckling  of  a  ring  under 
uniform  pressure  there  is  no  reason  to  use  a  finer  spacing  in  any  local  region. 
A  constant  mesh  spacing  gives  the  more  efficient  -  odel  and  the  results  for 
rings  with  the  mesh  indicated  in  Figure  12  are  shown  only  for  damns tration 
of  the  effect  of  variable  spacing. 
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Same  computed  buckling  loads  cure  shewn  in  Table  5.  Hie  second  order  accuracy 
is  maintained.  A  Richardson  extrapolation  based  on  the  results  for  9  and  12 
elements  over  the  90°  arch  leads  to  PCR  =  0.99997.  Comparison  with  the  results 
in  Table  4  indicates  that  the  results  are  somewhat  less  accurate  than  those 
obtained  with  the  finest  of  the  two  s pacings  used  and  somewhat  better  than 
those  obtained  with  the  coarser  spacing.  The  prebuckling  solution  for  the 
variable  mesh  with  6  elements  shows  a  rather  large  error  in  lateral  displacement 
although  the  error  in  the  buckling  loads  is  only  0.43%.  Hie  computed  displace¬ 
ments  vary  from  about  7%  above  to  about  7%  below  the  analytic  solution 
2 

(w  =  pR  /EA) ,  being  too  snail  in  the  area  with  coarse  spacing. 


Symmetry 


Figure  12.  Ring  with  Variable  Mesh 
TABLE  5 

BUCKLING  OF  RINGS  WITH  VARIABLE  MESH 

Number  of 

Elements,  i  D  Error 


Section  IX 


PLANE  STRESS  PROBLEMS 


In  plane  stress  problems  the  displacement  components  u,  v  are  functions  of 
two  spatial  coordinates  x,  y.  The  strain  energy  is  defined  in  terms  of  the 
first  order  derivatives  u,^,  u,y,  v,x,  v,  .  Consequently,  the  Ritz  theory  require 
that  the  functions  u  and  v  cure  continuous  over  the  domain.  The  simplest 
element  for  analysis  of  plane  stress  problems  is  the  triangular  element. 

A  Taylor  series  in  two  dimensions  including  derivatives  up  to  the  first  order 
contains  three  terms.  Consequently,  a  displacement  field  in  which  the  first 
order  derivatives  cure  of  first  order  accuracy  can  be  determined  from  three 
nodal  function  values.  Hie  obvious  choice  then  is  an  element  in  which  the 
degrees  of  freedom  are  represented  by  the  values  of  u  and  v  at  the  comers. 
Two  adjacent  elements  have  identical  values  of  the  displacement  components  at 
both  ends  of  the  interface  and  these  components  vary  linearly  with  the  space 
coordinates.  Thus,  the  values  of  u  and  v  in  the  two  elements  are  identical 
over  the  entire  interface,  i.e. ,  u  and  v  are  continuous  over  the  domain. 
Elements  that  satisfy  the  continuity  requirements  of  the  Ritz  theory  are 
referred  to  as  conforming  elements. 

With  linear  variation  of  the  displacements  all  strains  (in  linear  analysis)  are 
constant  within  an  element.  The  triangular  6-degree-of -freedom  elenent  is 
usually  referred  to  as  the  "Constant  Strain  Triangle"  or  the  CST-element.  As 
the  strain  is  constant  throughout  the  element  it  is  obviously  sufficient  to 
use  integration  with  one  Gaussian  point,  i.e.,  the  c.g.  of  the  triangle.  This 
point  is  also  the  stress  window  where  the  strain  is  expressed  with  a  second 
order  accuracy.  In  a  range  of  sufficiently  snail  elements  the  error  in  analysis 
with  CST  elements  varies  with  the  square  of  the  mesh  spacing. 

Equivalent  to  the  Taylor  series  approach  is  the  use  of  linear  shape  functions. 

Rar  triangular  elements  it  is  convenient  to  express  these  shapefunctions  in 
terms  of  so  called  area  coordinates  as  discussed,  for  example  in  Reference  12. 

The  definition  of  such  coordinates  is  illustrated  in  Figure  13. 


Figure  13.  Area  Coordinates  fear  Triangular  Elements 


The  location  of  a  point  P  inside  the  triangle  is  uniquely  determined  if  two 
of  the  areas  of  the  subtriangles  A^  are  known.  Consequently  these  areas,  or 
rather  the  areas  normalized  with  respect  to  the  total  area  A  of  the  elanent , 
can  be  used  as  space  coordinates. 


e  =  Aj/A  ,  n  =  A^/A  ,  c  =  Aj/A 


(48) 


As  the  relation  c  +  n  +  ?  =  1  holds,  the  Cartesian  coordinates  of  point  P 
can  be  written  in  the  form 


X1  x2 


3 

*3 


(49) 
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The  same  interpolation  may  be  used  for  the  displacement  components  u  and 
v,  i.e., 

:! 

where  u.v^  represent  the  discrete  values  of  the  displacement  components  at 
earner  i.  The  inversion  of  Equation  (49)  is  needed  for  expression  of  the 
location  of  Gaussian  integration  points  in  terms  of  the  area  coordinates. 
The  inverse  is  given  in  Reference  12. 


In  order  to  raise  the  order  of  accuracy  of  the  triangular  plane  stress 
element  it  is  necessary  to  add  another  three  freedoms  for  each  displacsnent 
component  (there  are  three  second  order  derivatives) .  A  natural  choice  then 
is  to  add  the  displacements  at  midside  nodes  on  each  of  the  element  boundaries. 
Figure  14  shows  a  triangular  element  with  midside  nodes  and  indicates  nodal 
values  of  the  area  coordinates. 
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Figure  14.  Linear  Strain  Triangle 
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The  relation  between  Cartesian  and  area  coordinates  is  given  by: 
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(51) 


Using  the  corresponding  interpolation  for  the  displacement  components  we 
obtain  a  simple  expansion,  equivalent  to  the  result  of  a  quadratic  Taylor 
series  approach. 


That  is 


(52) 


where  each  of  the  six  shape  functions 

Pi  =  U(1  -  20  ,  n(l  -  2n)  ,  s(  1  -  2n)  ,  4£n  ,  4nc  ,  4a) 
are  equal  to  unity  at  one  of  the  six  nodes  and  zero  at  the  others. 


Fran  Equations  (52)  and  (53)  it  follows  that  the  displacement  components 
are  quadratic  with  respect  to  the  spatial  coordinates.  Since  tvo  adjacent 
elements  have  three  freedoms  (for  each  displacement  component)  in  cannon 
the  linear  strain  triangle  (LST  element)  is  conforming  and  the  convergence 
requirements  of  the  Ritz  procedure  are  satisfied. 


(53) 


With  the  same  nodal  pattern  we  can  define  the  displacements  inside  the  element 
by  use  of  four  constant  strain  triangles.  Clearly  this  gives  a  less  accurate 
representation  of  the  strain  than  the  linear  variation  obtained  by  use  of  the 
shape  functions  in  Equation  ( 51) .  We  would  expect  therefore  that  the  higher 
order  representation  (the  linear  strain  triangle)  should  lead  to  a  more 
efficient  analysis. 


Triangular  plane  strain  elements  with  fourth  order  accuracy  require  at  least 

a  displacement  field  based  on  ten  degrees  of  freedom  (there  are  four  third 

order  derivatives) .  Such  elements  can  be  derived  (see  Reference  12,  for 

exanple)  by  introduction  of  a  node  internal  to  the  triangle  and  use  of  two 

midside  nodes  on  each  side.  Alternatively  freedoms  at  midside  nodes  can  be 

substituted  by  displacement  derivatives  at  the  earner  nodes.  This  seems  to  be 

undesirable  because  it  involves  use  of  normal  strains  u,  and  v,  as  freedoms. 

x  y 

These  are  generally  not  defined  on  boundaries. 


Three  sided  elements  with  curved  boundaries  can  be  defined  in  the  same  way  as 
the  triangular  elements.  Generally  such  elements  will  not  be  free  from  strain 
energy  due  to  a  rigid  body  displacement  (compare  the  discussion  of  curved 
beam  elements  above),  if  the  element  is  rigidly  displaced  in  its  plane  a  strain 
free  configuration  is  possible  only  if  the  shape  functions  used  for  displace¬ 
ment  can  exactly  represent  the  initial  shape  of  the  element  boundary.  In 
other  cases  an  error  is  introduced  which  disappears  with  diminishing  gridsize. 
The  order  of  the  error  depends  on  the  order  of  the  power  series  representing 
the  displacements.  One  way  to  avoid  self-straining  is  referred  to  as 
isoparametric  representation .  In  this  procedure  the  same  shape  functions  are 
employed  to  describe  (approximately)  the  element  boundaries  as  well  as  the 
inplane  displacement  configuration .  For  example,  with  1ST  elements  the 
element  boundaries  can  be  approximated  by  polynomials  up  to  the  second 
order. 


The  relative  efficiency  of  linear  strain  triangles  versus  constant  strain 
triangles  is  illustrated  in  Figure  15.  A  section  of  an  annular  plate  as  shewn 
in  the  figure  is  clamped  at  Side  1.  The  curved  sides  are  free  and  Side  3 
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Figure  15.  Deformation  of  Annular  Plate 


46 


is  subjected  to  a  displacement  of  0.1  in.  and  a  rotation  of  0.01  radians  as 

shewn.  The  figure  shews  the  error  in  the  total  reaction  force,  normal  to 

Side  1.  It  appears  that  for  the  same  accuracy  in  the  results  the  analysis 

with  constant  strain  triangles  leads  to  acmputer  cost  that  is  sane  twenty 

7 

times  higher  than  it  is  with  linear  strain  triangles.  With  E  =  10  psi, 
v  =  0.3  and  the  plate  thickness  0.1  in.  and  a  17  by  17  grid  the  total  reaction 
is  52905  lbs.  The  linear  strain  triangle  (or  quadrilateral)  should  only  be 
used  as  the  membrane  part  of  a  plate  element  and  its  use  should  be  restricted 
to  situations  in  which  accurate  representation  of  the  bending  behavior 
requires  a  grid  fine  enough  for  analysis  of  the  membrane  behavior  with 
constant  strain  elements. 

Quadrilateral  plane  stress  elotvents  can  conveniently  be  defined  by  combination 
of  two  or  more  triangular  elements.  In  case  of  a  linear  analysis  the  freedoms 
corresponding  to  nodes  that  became  internal  to  the  element  can  then  be 
eliminated  by  condensation.  That  is,  they  can  be  eliminated  through  energy 
minimization  on  the  element  level  and  do  not  appear  as  freedoms  in  the  final 
trial  functions. 

It  is,  of  course,  also  possible  to  derive  quadrilateral  elements  directly 
in  the  same  way  as  was  done  for  triangular  elements.  In  order  that  the 
displacements  at  nodes  be  compatible  with  that  of  adjacent  elements  it  is 
necessary  that  the  displacements  at  all  four  comers  be  included  as  degrees 
of  freedom.  Consequently,  it  is  generally  not  possible  to  restrict  the 
displacement  pattern  to  one  that  represents  constant  strain.  Using  the 
Taylor  series  approach  we  must  include  the  function  itself,  the  two  first 
order  derivatives  and  at  least  one  of  the  second  order  derivatives.  To  avoid 
directional  bias  we  must  choose  the  mixed  derivative  on  the  right  hand  side 
of  Equation  (4),  i.e.,  the  bilinear  expansion.  For  convenience  we  introduce 
the  coordinate  system  £,  n  with  the  property  that  £  equals  -1  and  +1  respec¬ 
tively  on  two  opposite  sides  and  tj  similarly  equals  -1  and  +1  on  the  other 
two  sides  as  shewn  in  Figure  16.  The  coordinates  £  and  n  are  expressed  in 
terms  of  the  x  and  y  coordinates  at  the  element  comers  by 


47 


where 


P1  =  J(1  -  a  (1  -n)  P3  =  j(l  +  £)  (1  -n) 

P2  =  £(1  +  £)  <1  -n)  P4  =  J(1  -  0  (1  -n) 

Similarly  a  displacement  component  can  be  defined  by 

4 

u  =  E  uipi  (55) 

i  =  1  11 

where  the  vu  are  the  nodal  values  of  the  displacanent  u.  In  this  case 
identical  interpolating  polynomials  are  used  for  displacements  and 
coordinates  (isoparametric  representation) .  The  mapping  (Figure  16) 
from  the  x,y  to  the  £,n  space  is  referred  to  as  isoparametric  mapping. 

Two  adjacent  element  have  common  displacements  at  the  endpoints  of  the 
interface.  Since  either  £  or  n  is  constant  along  any  element  boundary 
the  displacement  field  is  linear  along  boundaries  and  an  element  based 
on  these  shape  functions  is  conforming. 

In  order  to  achieve  second  order  accuracy  in  the  definition  of  the  first  order 
derivatives  anywhere  inside  the  element  we  must  include  at  least  six  degrees 
of  freedom.  Directional  bias  will  result  unless  the  number  of  nodes  is 
divisible  by  four,  a  possible  node  at  the  element  midpoint  not  counted. 
Consequently  the  higher  order  element  should  have  8  (or  9)  nodes.  It 
seems  logical  to  add  a  node  at  the  middle  of  each  side. 

The  higher  order  representation  of  the  displacements  allows  parabolic 
approximation  of  curved  element  boundaries  without  self- straining. 


In  general,  third  order  accuracy  in  the  first  order  derivatives  requires 
that  at  least  ten  degrees  of  freedom  are  utilized  in  the  definition  of  the 
displacement  field.  However,  in  the  one-dimensional  case,  the  three  point 
scheme  for  inplane  displacements  shown  in  Figure  8  results  in  an  expression 
for  the  first  order  derivative  with  a  third  order  error  at  the  two  Gaussian 
point.  Therefore,  we  can  expect  a  third  order  accuracy  in  the  element  with 
nine  nodes.  These  linear  strain  quadrilateral  elements  are  referred  to  as 
the  1SQ8  and  LSQ 9  elements.  With  two  inplane  displacement  components  these 
elements  have  16  and  18  degrees  of  freedoms,  respectively. 
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Section  x 

PLATE  AND  SHELL  ELEMENTS 

The  energy  functional  based  on  plate  or  shell  theory  contains  changes  of  curva¬ 
ture  of  the  reference  surface,  i.e.  second  order  derivatives  of  the  lateral 
displacements.  Consequently,  the  trial  functions  as  well  as  their  first  order 
spatial  derivatives  must  be  continuous  across  element  boundaries  in  order  that 
convergence  with  decreasing  grid  size  to  the  correct  solution  will  follow 
from  the  Ritz  theory.  Efforts  to  define  efficient  and  conforming  plate  or 
shell  elements  has  led  to  a  great  proliferation  of  bending  element  configura¬ 
tions. 

In  linear  analysis  of  flat  homogenous  plates,  with  the  midsurface  as  reference 
surface,  membrane  and  bending  behavior  cure  uncoupled .  The  bending  energy  in 
the  plate  depends  exclusively  on  the  lateral  displacement  w  and  the  membrane 
energy  exclusively  on  the  inplane  displacements  u  and  v.  in  that  case  it  is 
possible  to  superimpose  bending  elements  onto  the  membrane  (plane  stress) 
elements  discussed  in  the  preceding  paragraphs.  In  nonlinear  analysis  the 
membrane  energy,  due  to  stretching  of  the  middle  surface,  depends  on  the 
lateral  displacement  pattern  and  in  analysis  of  curved  elements  the  bending 
energy  is  a  function  of  inplane  as  well  as  lateral  displacements. 

Bending  and  membrane  elements  then  have  ocmton  degrees  of  freedom  but  as  long 
as  the  constitutive  equations  for  the  shell  wall  do  not  introduce  manbrane- 
bending  coupling  it  is  still  possible  to  superinpose  independent  balding  and 
membrane  elements.  Most  problems  involved  in  the  development  of  finite  element 
configurations  for  shell  or  plate  analysis  are  independent  of  this  coupling. 
Therefore,  it  is  helpful  to  discuss  same  aspects  of  these  configurations  in  the 
less  ccmple  framework  of  pure  bending  analysis. 

Plate  bending  elements  generally  include  as  degrees  of  freedom  rotations 
(displacement  derivatives)  as  well  as  lateral  displacements .  in  a  triangular 
element  the  number  of  degrees  of  freedom,  must  be  divisible  by  three,  those 
at  a  possible  node  at  the  midpoint  of  the  element  not  counted,  as  otherwise 
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directional  bias  would  occur.  Second  order-derivatives  must  be  determined 
at  least  to  a  first  order  accuracy.  As  a  consequence  the  lowest  possible 
number  of  degrees  of  freedom  is  six. 


The  element  in  Figure  17  shewn  in  a  Cartesian  system  x,  y,  z  then  represents 
the  simplest  possible  triangular  bending  element.  Use  of  a  Taylor  series 
approach  leads  to  a  uniquely  determined  complete  quadratic  representation  of 
the  lateral  displacement  field,  that  is 


2  2 
w  =  a  +  bjX  +  b2y  +  c^x  +  c^xy  +  C22^ 


W3 


(56) 


The  changes  of  curvatures  as  second  order  spatial  derivaties  of  w  are  all 
constant  over  the  element  surface.  Along  a  side,  between  nodes  1  and  2 
(y  =  0)  for  example,  the  rotation  around  the  element  boundary  is  of  the  form 

w,y  =  bl  +  c12* 

That  is,  the  slope  varies  linearly.  However,  two  adjacent  element  have 
only  one  rotational  freedom  in  common.  Therefore,  the  constant  bending 
element  is  nonconforming.  Similar  problems  are  encountered  when  triangular 
bending  elements  of  higher  order  are  developed.  Elements  with  nine  and 
twelve  freedoms  are  shewn  in  Figure  18. 


Figure  18.  Triangular  Elements  with  Nine  (a)  and  Twelve  (b) 

Degrees  of  Freedom 

A  complete  cubic  can  be  determined  by  use  of  ten  degrees  of  freedom.  Thus  the 

displacement  field  in  the  nine  degrees  of  freedom  element,  can  be  determined 

by  use  of  a  cubic  Taylor  series  only  after  a  constraint  is  introduced.  For 

2  2 

example,  the  coefficients  for  the  x  y  or  xy  terms  may  be  required  to  be 
identical.  In  that  case  the  expression  for  the  rotation  around  the  element 
boundary  contains  terms  which  are  second  order  in  a  coordinate  along  the  boundary. 
As  only  two  rotational  freedoms  sure  oonmon  between  adjacent  elements  rotational 
nonconformity  is  allowed. 

With  twelve  degree  of  freedom  element  two  quartic  terms  are  included.  Either 
4  4  3  3 

of  the  pairs  x  ,  y  or  x  y,  xy  can  be  included  but  in  any  case  the  edge 
rotation  becomes  cubic  in  terms  of  the  spatial  coordinate  along  the  boundary. 

As  four  coefficients  are  required  to  determine  uniquely  a  cubic  and  adjacent 
elements  only  have  three  rotational  freedoms  in  common,  the  element  is  non¬ 
conforming. 


It  is  possible  to  eliminate  slope  nonconf  ormi ty  through  introduction  of  displace¬ 
ment  constraints  or  to  minimize  its  effect  through  introduction  of  penalty 
functions.  Hie  latter  approach  can  be  applied  to  the  twelve  degrees  of  freedom 
element.  A  fictitious  term 

AU  =  Cx  a^  +  C2  a^3  (58) 

is  added  to  the  energy,  where  a^  and  a-^  are  the  coefficients  in  the  x3y 
and  xy^  terms  while  and  C 2  are  constants  chosen  large  enough  to  provide 
a  sufficient  penalty  on  the  strain  energy  due  to  the  quartic  terms. 

Introduction  of  constraints  (or  penalty  functions)  to  achieve  conformity 
lowers  the  order  of  the  accuracy  of  the  local  approximation.  Since  the  six 
degree  of  freedom  bending  triangle  is  of  first  order  accuracy  only  it  cannot 
be  made  conforming  by  introduction  of  constraints. 

A  conforming  triangular  bending  element,  introduced  by  Clough  and  Tocher  in 
Reference  13,  is  based  on  partition  into  three  subtriangles.  This  elanent  is 
discussed  in  sane  detail  in  Reference  12.  Each  subtriangle  has  10  degrees 
of  freedom  as  shown  in  Figure  19  allowing  a  complete  cubic  representation. 

The  directional  bias  in  the  subtriangles  is  eliminated  when  they  are  combined 
into  one  element.  Since  the  internal  sub  triangle  boundaries  do  not  include 
midpoint  rotations  as  freedoms  it  is  necessary  for  rotational  compatibility 
to  constrain  the  element  so  that  the  rotation  varies  linearly  along  these 
boundaries.  The  combined  element  has  3  freedoms  at  each  corner,  3  freedoms 
at  the  element  midpoint  and  3  midside  rotational  freedoms  for  a  total  of 
15  freedoms.  For  linear  analysis  the  freedoms  at  the  element  midpoint  can  be 
eliminated  on  the  element  level  through  condensation.  While  the  introduction 
of  contraints  reduces  the  order  of  accuracy,  condensation  does  not  have  any 
effect  on  the  final  solution. 


Figure  19.  The  Clough-Tbcher  Triangular  Bending  Element 


Section  XI 


QUADRILATERAL  BENDING  ELEMENTS 

The  number  of  freedoms  in  a  quadrilateral  element  must  be  divisible  by  four 
in  order  that  directional  bias  be  avoided  (possible  freedoms  at  an  elanent 
midpoint  node  discarded) .  A  bending  element  mist  as  a  minimum  include  all  six 
terms  up  to  the  second  order.  Consequently  the  simplest  possible  quadrilateral 
element  would  have  8  degrees  of  freedoms.  Experience  with  constant  strain 
elements  for  plane  stress  analysis  indicates  that  such  an  element  is  not 
likely  to  be  very  efficient. 

The  first  element  for  analysis  of  the  bending  of  flat  plates  was  a  rectangular 

element  with  the  freedom  pattern  illustrated  in  Figure  20.  As  the  element  has 

12  degrees  of  freedom,  the  displacement  pattern  can,  as  in  the  triangular  element 

with  the  same  number  of  freedoms  be  represented  by  a  complete  cubic  and  two 

3  3 

quartic  terms,  for  example,  x  y  and  yx  .  Such  an  elanent  was  introduced 
(Reference  14)  and  successfully  applied  in  plate  bending  analysis  before  the 
finite  element  method  was  established  as  a  form  of  the  Ritz  procedure.  We 
refer  to  this  element  with  the  notation  QB12  (quadrilateral  bending,  12 
freedoms) .  Shape  functions  for  this  element  equivalent  to  the  Taylor  series 
solution  are  presented  in  Reference  15. 


After  the  finite  element  method  was  given  a  solid  theoretical  foundation,  the 
wisdom  in  the  use  of  nonconforming  element  was  questioned.  Efforts  were  made 
to  improve  the  performance  of  the  bending  elements  through  introduction  of 
displacement  constraints.  As  the  order  of  accuracy  of  the  local  approxima¬ 
tions  thus  was  reduced  these  efforts  were  not  particularly  successful. 

Generally  the  nonconforming  elements  led  to  better  results  within  the  range 
of  engineering  accuracy. 

In  the  pursuit  of  more  efficient  conforming  bending  elanents  Clough  and  Felippa 
derived  quadrilateral  elements  through  the  combination  of  the  Clough-Tocher 
triangular  elements.  Either  two  or  four  triangles  are  combined  to  form 
quadrilateral  elements  and  for  linear  analysis  internal  freedoms  are  eliminated 
through  condensation.  After  condensation  the  elements  have  12  degrees  of 
freedom  and  we  refer  to  them  as  the  CF12,2  or  CF12,4  elanents  depending  on 
the  number  of  sub triangles  vised.  With  three  rotational  freedoms  along  each 
side  and  only  quadratic  representation  of  the  rotations (cubic  in  the  displace¬ 
ments)  the  element  is  conforming.  The  disadvantage  with  the  element  is  that 
the  subtriangle  approach  makes  it  necessary  to  use  a  large  number  of  integration 
points  resulting  in  excessive  ccrrpufcer  time. 

A  conforming  bending  element  was  presented  by  Bogner,  Fox  and  Schmit  in 

Reference  16.  A  theory  for  curved  shells  is  utilized  in  the  derivation  of  the 

stiffness  matrix.  The  strain  energy  due  to  rigid  body  displacement  is  exactly 

zero  (independent  of  grid  size)  only  for  flat  rectangular  elements.  Shape 

functions  are  used  for  interpolation  to  ensure  conformity.  The  only  version 

of  the  element  that  can  be  seriously  considered  is  one  BFS16  with  16  degrees 

of  freedom  in  which  the  twist,  w,  ,  has  been  added  as  a  freedom  at  each 

xy 

node.  While  the  interpolating  polynomials  contain  terms  up  to  the  sixth  order 
the  rotation  along  the  shell  boundary  is  cubic. 

Seme  results  obtained  with  this  element  are  presented  in  Reference  16.  A  flat 
plate  20  x  20  in  0.1  in.  thick  with  clamped  edges  is  loaded  by  a  uniform 
normal  pressure  of  0.2  psi.  The  material  is  characterized  by  E  =  107  psi, 


v  =  0.3.  The  results  are  shewn  in  Table  6  together  with  corresponding  results 
based  on  analysis  with  the  CF12,2  element. 


Table  6 

BENDING  OF  SQUARE  PLATE 

Number  of  Elements  Midpoint  Displacement  (in) 

(on  one  quarter  of  the  plate)  CF12,2  BFS  16 


2  x  2 

.035592 

.040475 

3x3 

.039713 

.040482 

4  x  4 

.040274 

.040487 

5  x  5 

.040409 

6x6 

.040453 

10  x  10 

.040485 

The  convergence  with  the  BFS  16  element  is  exceptionally  good.  For  example, 
the  solution  based  on  16  such  elements  is  as  good  as  one  based  on  100 
CF12 , 2  elements .  Still  it  appears  that  the  BSF16  element  has  seen  little 
use.  One  reason  may  be  that  it  is  self -straining  unless  flat  and  rectangular. 
A  more  important  disadvantage  is  that  the  twist  w^  is  used  as  a  degree  of 
freedom.  This  displacement  parameter  cannot  readily  be  defined  on  shell 
boundaries. 


Section  XII 


HYBRID  METHODS 

The  finite  element  method  was  originally  developed  as  a  direct  extension  of 
the  methods  for  analysis  of  statically  determined  structures.  In  these 
methods  internal  forces  were  introduced  as  unknowns  rather  than  displacements 
(Reference  17,  for  example) .  The  derivation  of  equilibrium  equations  was 
based  on  the  principle  of  minimum  complementary  energy.  While  the  displace¬ 
ment  method  has  became  the  dominant  procedure  in  finite  element  analysis  the 
so  called  force  method  still  shews  sane  signs  of  life.  A  plate  element  based 
on  the  principle  was  recently  presented  in  Reference  18. 

The  difficulties  encountered  in  satisfying  continuity  conditions  at  boundaries 
between  adjacent  elements  for  shell  or  plate  bending  elements  has  motivated 
the  pursuit  of  other  avenues,  such  as  the  development  of  mixed  or  hybrid 
methods  in  which  both  displacement  and  force  fields  are  assumed.  The  hybrid 
method  developed  by  Pian  and  Tong  (Reference  19)  represents  an  important 
contribution  in  finite  element  technology.  In  this  method  the  displacement 
freedoms  are  not  directly  used  for  definition  of  the  displacements  inside  the 
element.  A  stress  field  with  determined  coefficients  is  assumed  in  addition 
to  the  nodal  displacements  on  the  element  boundary.  The  unknown  coefficients 
in  the  stress  field  are  determined  through  energy  minimization  on  the  element 
level  with  the  nodal  displacements  defining  the  boundary  conditions.  The 
element  stiffness  matrix  so  obtained  is  based  on  the  nodal  displacement 
freedoms. 

As  the  principle  of  canplementary  strain  energy  is  to  be  applied  the  strain 
energy  is  written  in  berms  of  stresses,  i.e. , 

U  =  j  f  { a )T  [N]  {a}  dv  (59) 

v 

where  [N]  is  the  matrix  of  compliance.  That  is  the  inverse  of  the  stiffness 
matrix. 

{a}  =  [C]  U)  ,  N  =  [C I-1  (60) 
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The  stress  vector  is  defined  in  terms  or  a  number  of  undetermined  coefficients 
(6)  (stress  freedoms) 


{a}  =  [P]  {&}  (61) 

where  the  elements  of  the  matrix  IP]  are  functions  (polynomials)  of  the 
spatial  coordinates  similar  to  the  displacement  shape  functions  discussed 
above.  The  strain  energy  in  terms  of  stresses  then  can  be  written  in 
the  form 

U  =  |  (6}T  [H]  {0}  (62) 

where 

H  =  /  [P]T  [N]  (P]  dV  (63) 

V 

The  displacements  along  element  boundaries  are  expressed  in  terms  of  the 
nodal  displacement  freedoms.  Forces  along  the  boundaries  are  obtained 
by  substitution  of  the  appropriate  values  of  the  spatial  coordinates 
in  the  expression  for  the  stress  field, Equation  (61) .  After  integration 
along  element  boundaries  the  work  done  by  these  boundary  forces  is 
readily  obtained 


fi  =  (0}T  [T]  {q}  (64) 

where  T{q}  represents  an  interpolation  of  the  displacements  along 
element  boundaries  in  terms  of  the  vector  of  nodal  freedoms  {q>. 

Minimization  of  the  conplementary  energy 


V  =  U  -  n 


(65) 


60 


leads  to 


(6)  =  [H]"1  IT]  {q}  (66) 

and  the  total  strain  energy  can  be  written  as 

U  =  \  {q}T  [K]  {q}  (67) 


where 


K  =  [T]T  [H]"1  [T]  (68) 

is  the  element  stiffness  matrix  corresponding  to  the  displacement  freedoms  q. 

While  the  use  of  forces  as  freedoms  leads  to  an  underestimate  of  the  strain 
energy,  the  displacement  method,  provided  the  rules  of  the  Ritz  theory  are 
followed,  overestimates  the  strain  energy.  The  hybrid  method  as  presented 
above  gives  results  that  are  bracketed  by  those  obtained  from  force  and 
displacement  methods.  Consequently,  the  convergence  in  many  cases  is 
very  good.  A  rectangular  bending  element  of  this  type  was  presented  in 
Reference  19.  This  element  has  12  displacement  freedoms  and  up  to  23 
coefficients  in  the  stress  field.  A  general  quadrilateral  version  of  the 
element  is  used  in  the  SPAR  computer  program  (Reference  20) .  We  refer  to 
this  element  by  the  notation  QH12  (quadrilateral  hybrid,  12  freedoms) . 
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Section  XIII 


USE  OF  NOTCONFOFMING  ELEMENTS 

After  the  finite  element  analysis  was  given  a  firm  theoretical  foundation  as  a 
form  of  the  Ritz  procedure,  the  question  of  possible  convergence  of  nonconform¬ 
ing  elements  has  been  a  subject  of  considerable  interest  (References  6  and  21, 
for  example) .  The  prooess  of  integration  of  the  energy  over  the  individual 
elements  followed  by  a  sirmation  over  all  elanents  leaves  out  possible  work 
done  by  forces  on  the  element  boundaries.  For  conforming  elements  this  work 
vanishes  because  contributions  from  adjacent  elements  cancel  one  another. 

If  the  elements  are  nonconforming  the  work  corresponding  to  the  forces  on  the 
boundary  and  the  discontinuity  in  displacement  (or  rotation)  is  left  out  of  the 
energy  balance.  The  procedure  can  converge  to  the  correct  solution  only  if 
this  work  vanishes  with  diminishing  grid  size. 

In  Reference  22  Irons  and  his  coworkers  suggested  a  simple  test  to  be  applied 
to  nonconforming  elements.  This  test  has  later  become  known  as  the  patch  test 
and  it  appears  that  the  passing  of  this  test  is  a  sufficient  but  not  always 
necessary  convergence  requirement.  As  first  presented  the  patch  test  required 
that  for  ary  patch  of  elements  subjected  to  displacements  corresponding  to 
constant  strain  (constant  curvature  for  bending  elements)  on  its  boundary,  the 
solution  everywhere  in  the  interior  should  be  accurate,  i.e.,  the  strain  is 
constant  throughout  the  patch. 

In  Reference  6  Strang  and  Fix  show  that  the  patch  test  is  equivalent  to  the 
requirement  that  the  work  corresponding  to  displacement  discontinuities  along 
element  boundaries  vanishes.  A  rectangular  element  developed  by  Wilson 
(Reference  23)  is  used  for  demonstration.  Wilson  raises  the  order  of  the 
hi) inear  plane  stress  elenents  (Figure  16)  by  addition  of  the  two  so  called 
bubble  modes 


P5  =  (1  -  £2)  and  Pg 


(1  -  n2) 


(69) 
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These  are  clearly  nonconforming  nodes,  independent  on  the  nodal  displacement 
components.  The  formulation  was  suggested  for  rectangular  elements  only  and 
for  such  elements  it  is  shewn  in  Reference  6  that  the  addition  of  the  noncon¬ 
forming  bubble  nodes  results  in  a  considerable  improvement  in  the  convergence 
rate. 

It  is  shown  in  Reference  6  that  "the  energy  of  the  discontinuity"  vanishes 
identically  vfoen  the  integration  is  carried  over  the  entire  boundary.  For  a 
general  quadrilateral  element  on  the  other  hand,  the  integral  over  the  boundary 
can  be  represented  as  a  constant  times  the  difference  in  the  lengths  of  opposite 
sides. 

Consequently,  the  contribution  to  the  energy  from  one  element  is  proportional  to 

the  grid  spacing.  Hcwever,  the  total  length  of  all  element  boundaries  is  pro- 
2 

porticnal  to  1/h  .  Therefore,  it  does  not  seem  likely  in  any  given  case  that 
analysis  based  on  ncnrectangular  elements  with  bubble  modes  will  converge  to 
the  correct  solution.  A  modification  of  this  approach  presented  in  Reference  24 
makes  it  possible  to  use  the  bubble  inodes  also  with  ncnrectangular  elements. 

A  relatively  simple  and  often  used  nonconforming  bending  element  is  the  quadri¬ 
lateral  12  degree  of  freedom  element  QB12  (see  Figure  20) .  For  this  element  the 
lateral  displacement  component  can  be  expressed  in  the  form 

"  -  Co  +  Clx  *  C2*  +  C3x2  +  *  C5*2  (70) 

+  Cgx3  +  G^y  +  Cgxy2  +  C9y3 

+  C10x3y  +  Cnxy3 


The  element  edge  rotation  is  represented  by  the  normal  derivative,  that  is,  if  the 
side  makes  an  angle  a  with  the  y-axis 

3w/3n  =  3w/3x  cos  a  +  3w/3y  sin  a 


1 Vro  adjacent  elements  have  matching  slopes  at  the  end  points  of  the  ocmcn 
boundary.  Consequently  terms  of  lower  than  second  order  in  the  coordinate 
along  the  boundary  do  not  cause  any  slope  discontinuity.  The  mismatch  in 


slope  between  adjacent  elements  is  of  the  form 


A?  =  C£2  +  higher  order  terms  (72) 

To  obtain  the  energy  of  the  discontinuity  we  multiply  the  slope  discontinuity 
by  the  average  bending  moment  Mq  (higher  order  terms  are  irrelevant)  and 
integrate  the  product  over  the  length  h  of  the  interface.  The  result  is 
of  the  form  CMQh^.  As  the  total  length  of  element  boundaries  is  proportional 
to  1/h2,  the  total  energy  of  the  discontinuity  is  proportional  to  the  grid 
spacing.  We  expect  therefore  that  the  use  of  the  QB12  element  will  lead  to 
convergence  of  first  order,  i.e.,  E  =  0(h). 

If  the  element  is  rectangular  it  passes  the  patch  test  as  posed  in 
Reference  6.  Consequently,  any  displacement  pattern  with  constant  curvature 
is  exactly  represented,  independently  of  grid  size,  i.e.,  the  aocuracy  is 
at  least  of  second  order.  Ihe  general  quadrilateral  element  passes  a  milder 
form  of  the  test,  requiring  that  the  error  in  energy  introduced  by  the 
nonconforming  inodes  is  proportional  to  same  pcwer  n>0  of  the  gridsize. 

The  performance  of  the  QB12  element  is  illustrated  here  in  a  few  examples. 

A  simply  supported  square  plate  5x5  in.  and  0.1  in.  thick  is  made  of  a 

7 

material  with  E  =  10  psi  and  v  =  0.3.  The  plate  is  subjected  to  unidirec¬ 
tional  compression.  In  that  case  (see  Reference  25  for  example)  the 
critical  load  is 

2_3 

N_  =  4  -  V™-  7  =  1446.0958  lb/in.  (73) 

^  12(1  -  v)b2 

The  buckling  mode  may  be  assnned  symmetric  about  two  planes.  Hence,  in  the 
finite  element  analysis  it  is  sufficient  to  consider  a  quarter  model  of  the 
plate.  Some  results  are  shown  in  Table  7,  where  N  represents  the  number  of 
elements,  in  each  direction,  on  this  quarter  model. 
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TABLE  7 


BUCKLING  OF  SQUARE  PLATE 


No.  of  Elements 

QB12 

QJ2 

QH12 

EN4 

N  x  N 

Critical  Load  Error 

E 

Critical  Load  Error 

E 

2x2 

1391.19 

54.91 

220 

1451.04 

4.94 

79 

4x4 

1431.44 

14.66 

235 

1446.41 

.31 

79 

6x6 

1439.50 

6.60 

238 

1446.16 

.06 

78 

8x8 

1442.37 

3.73 

239 

10  x  10 

1443.71 

2.39 

239 

Richardson's  extrapolation  based  on  the  results  for  N  =  8  and  N  =  10  with  the 
assumption  of  second  order  convergence  leads  to  =  1446 . 100 ,  i.e. ,  to  an 
error  of  0.0007.8  percent.  It  can  be  little  doubt  then  that  the  use  of  this 
nonconforming  element  leads  to  a  second  order  convergence.  For  a  one  percent 
error  we  need  four  elements  per  quarter  wave.  Four  elements  per  half  wave 
gives  a  barely  acceptable  3.8  percent  error. 

Table  7  shows  also  results  obtained  with  the  hybrid  element  QK12.  In  this  case 
convergence  is  from  above  and  clearly  a  fourth  order  accuracy  is  obtained.  With 
four  elements  per  halfwave  the  error  is  as  snail  as  0.34  percent.  For  the  same 
grid  the  runtime  with  the  QS12  element  is  only  slightly  larger  than  with  QB12. 

The  problem  of  buckling  of  a  spherical  shell  under  uniform  external  pressure 
was  considered  in  a  study  of  the  convergence  behavior  of  the  nonrectangular 
QB12  element.  In  order  that  effects  of  approximations  in  the  inplane  displace¬ 
ment  field  be  eliminated,  such  displacements  were  suppressed  in  the  buckling 
mode.  The  computed  buckling  loads,  then  are  physically  meaningless  but  this  does 
not  mar  the  conclusions  about  the  numerical  behavior  of  the  bending  elemait.  A 
spherical  segnent  with  radius  10  in. ,  thickness  0.1  in. ,  Young's  modulus  10 ^  psi 
and  Poisson's  ratio  0.3  was  analyzed.  Circumferentially  the  segment  covers  9 
degrees.  In  the  direction  of  the  meridian  the  segnent  covers  the  range  36  to 
45  degrees  measured  from  the  apex.  This  segnent  buckles  in  a  mode  with  one  half 
wave  in  the  meridional  and  one  quarter  wave  in  the  circumferential  direction. 
Therefore,  it  was  assumed  that  the  buckling  mode  was  antisyrrmetric  with  respect 
to  one  of  the  sides  along  a  meridian.  Symnetry  conditions  were  applied  at  the 
other  three  sides.  Buckling  loads  are  shown  in  Table  8. 
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TABLE  8 


BUCKLING  OF  SPHERICAL  SHELL  SEGMENT 


Number  of 
Elements 

N  x  M 

Critical 

Load 

Error 

% 

4  x 

2 

9034.1 

12.82 

6  x 

3 

9704.4 

6.35 

8  x 

4 

9965.9 

3.83 

10  x 

5 

10091.9 

2.61 

12  x 

6 

10161.7 

1.94 

16  x 

8 

10232.1 

1.26 

20  x 

10 

10265.0 

0.95 

24  x 

12 

10283.3 

0.77 

28  x 

14 

10295.0 

0.65 

32  x 

16 

10303.5 

0.57 

It  appears  that  the  convergence  is  of  first  order  and  for  grids  as  fine  as 
24  x  12  or  better  the  first  order  term  dominates  the  error.  In  that  ^.ase  the 
buckling  load  would  be  10362.6  psi.  The  errors  given  in  the  table  are  based 
on  this  value.  The  convergence  in  this  case  is  much  worse  than  it  is  with 
rectangular  elements.  Eight  elements  per  halfwave  correspond  to  an  error  of 
3.8%,  and  for  a  one  percent  accuracy  we  need  20  elements  per  half  wave. 

The  results  from  two  linear  plate  bending  analyses  are  shown  in  Figure  21. 

In  one  case  a  quadratic  plate  with  clamped  edges  is  subjected  to  uniform 
lateral  pressure.  The  error  in  the  midpoint  displacement  is  shown  as  a  function 
of  the  number  of  elements  along  each  side  on  a  quarter  model  of  the  plate. 
Despite  the  nonconformity  the  displacement  is  underestimated. 

The  other  case  is  one  with  nonrectangular  elements.  The  annular  plate  is 
clamped  at  one  of  the  straight  edges  and  the  other  three  edges  are  free.  The 
plate  is  subjected  to  a  uniform  pressure.  The  error  in  lateral  displacement 
at  point  A  is  shown  in  the  figure.  As  a  result  of  rotational  nonconformity, 
convergence  is  now  from  the  opposite  direction,  i.e. ,  the  displacements  are 
overestimated.  However,  the  convergence  is  still  relatively  good 
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Figure  21.  Error  Versus  Gridsize 


presumably  because  inaccuracies  due  to  the  nonconformity  and  approxima^  ;on 
of  the  deformation  mode  tend  to  cancel  one  another. 


For  some  reason,  presently  unexplained,  the  bending  analysis  of  the  annular 
plate  leads  to  an  illconditioned  equation  system.  With  grids  13  x  13  or  finer 
it  was  necessary  to  apply  a  linear  refinement  procedure  in  order  to  obtain 
accurate  results  (on  CDC  175) .  This  problem  did  not  occur  in  the  buckling 
analysis  of  the  spherical  shell  segment.  However,  with  slew  convergence  for 
the  spherical  shell  and  conditioning  problems  in  the  annular  plate  analysis 
it  seems  advisable  to  avoid  use  of  9312  elements  that  deviate  substantially 
from  a  rectangular  plan  form.  The  rectangular  version  of  the  element  may 
often  be  unsuitable  for  buckling  analysis  because  the  critical  load  converges 
frcm  below.  This  restricts  the  possibilities  to  increase  the  grid  spacing 
in  areas  with  light  loading. 
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Section  XIV 


AHMAD-TYPE  ELEMENTS 

Three-dimensional  elements  can  easily  be  derived  by  direct  extension  of  the 
isoparametric  representation  far  the  plane  stress  case  discussed  above 
(see  Figure  15) .  Second  order  derivatives  (curvature  changes)  do  not  appear 
in  the  functional,  and  therefore  it  is  sufficient  for  convergence  that  the 
displacement  components  themselves  are  conforming.  There  is  no  requirement 
of  interelement  compatibility  of  displacement  derivatives  (rotations) .  In 
view  of  the  difficulties  encountered  in  the  derivation  of  elements  based  on 
the  theory  of  plates  and  shells,  it  appears  worthwhile  to  try  to  represent 
the  shell  by  a  single  layer  of  such  brick  elements  (Figure  22) . 


Figure  22.  Isoparametric  Brick  Element  with  20  Nodes 

Primarily  due  to  the  problem  with  illconditioned  matrices  this  approach  did 
at  first  not  meet  with  much  success.  The  thinner  the  shell,  the  less  suitable 
cure  the  brick  elements.  This  problan  was  circumvented  in  a  publication  by 


Ahmad,  Zienkiewicz,  and  irons  (Reference  26)  by  introduction  of  the  basic 
assumptions  of  the  shell  theory  into  the  formulation  of  the  brick  elenent . 

First  it  is  assumed  that  the  energy  due  to  stresses  in  the  direction  of  the 
shell  normal  can  be  disregarded.  The  assumption  that  normals  to  the  shell 
remain  straight  is  also  invoked.  As  a  consequence  it  is  possible  to  express 
the  displacement  components  at  nodes  on  the  shell  surface  in  terms  of  displace¬ 
ments  and  rotation  components  at  the  shell  reference  surface  by  use  of  simple 
transformation  matrices.  The  set  of  freedoms  referred  to  in  the  interface  with 
a  program  user  is  shown  in  Figure  23a  and  the  equivalent  set  of  basic  freedoms 
on  which  the  energy  expression  is  based  is  shown  in  Figure  23b. 
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a)  "User  Freedom*"  fc)  Basic  Freedoms 

Figure  23.  Freedom  Patterns  in  an  Ahmad  Type  Element 

The  user  freedoms  are  three  displacements  and  two  rotations  at  each  of  eight 
reference  surface  nodes  for  a  total  of  40  freedoms.  The  basic  freedom  pattern 
contains  two  freedoms  at  each  of  16  shell  surface  nodes  and  one  freedom  at 
each  of  8  midsurface  nodes,  again  for  a  total  of  40  freedoms. 


In  a  first  carder  shell  theory  it  is  assumed  that  normals  ranain  straight  and 
normal  to  the  reference  surface  in  the  deformed  configuration.  This  means 
that  the  theory  neglects  deforma ttions  due  to  transverse  shear  stresses. 

A  second  order  theory  approximately  includes  transverse  shear  effects  as  the 
normal  is  assumed  to  remain  straight  but  is  allowed  to  rotate  with  respect  to 
the  shell  reference  surface.  Hie  shell  assumptions  introduced  in  the  Ahmad 
element  correspond  to  those  of  a  second  order  theory. 

A  first  order  shell  theory  also  contains  the  assumption  that  the  ratio  of  thick¬ 
ness  to  shell  radius  is  snail  so  that  it  can  be  neglected  in  comparison  to  one. 
This  assumption  is  not  introduced  in  an  Ahmad  tvpe  element.  A  more  accurate 
description  of  the  shell  bending  deformation  is  therefore  obtained  and  together 
with  the  inclusion  of  transverse  shear  deformation  this  extends  the  applicabil¬ 
ity  of  an  analysis  with  Ahmad  type  elements  to  the  range  of  "moderately  thick" 
shells . 

The  inclusion  of  transverse  shear  deformation  introduces  a  special  problem. 

The  description  of  the  deformation  pattern  must  be  general  enough  to  inhibit 
the  development  of  excessive  transverse  shear  strain.  The  presence  of  such 
strains  leads  to  unsatisfactory  convergence  properties  with  the  early  Ahmad 
type  elements.  In  Reference  27  Pawsey  circumvents  this  problem  through 
introduction  of  reduced  integration.  Following  Reference  27  we  illustrate  the 
problan  in  the  one  dimensional  case.  We  consider  the  case  in  which  the  bend¬ 
ing  moment  varies  linearly  and  vanishes  at  the  midpoint  of  the  beam  element  as 
illustrated  in  Figure  24a. 


Figure  24.  Ahmad  Type  Beam  Elanent  with  Linearly  Varying  Bending  Moment 

The  correct  solution  with  the  rotations  +0,  -6,  +0  at  the  three  nodes  is  shown 
in  Figure  24b.  Since  the  shape  functions  are  not  sufficiently  general  to  allow 
this  deformation  pattern,  the  solution  represented  in  Figure  24c  is  instead 
obtained,  in  this  deformation  pattern  the  transverse  shear  contributes  a  very 
large  part  to  the  strain  energy.  However,  Pawsey  shews  that  the  correct  value 
of  the  shear  strain  is  obtained  at  the  Gaussian  points  with  two-point  integration. 
A  similar  problem  is  encountered  far  curved  beam  elements  subjected  to  a 
constant  bending  moment.  For  this  case  Pawsey  shows  that  a  spurious  normal 
membrane  strain  develops,  but  also  that  this  strain  vanishes  at  the  two 
Gaussian  points. 

It  was  established  thus  in  Reference  27  that  for  the  Ahmad  type  elanent, 
convergence  with  gridsize  is  considerably  better  if  we  use  a  reduced  integra¬ 
tion  schene.  For  the  plate  bending  element  this  would  mean  that  eight  (two 
sets  of  four  points)  integration  points  are  used  rather  than  the  eighteen 
points  that  would  be  required  for  accurate  integration  of  the  linear  part 
of  the  strain  energy.  As  noted  above  the  use  of  reduced  integration  can 


allow  the  development  of  strain-free  deformation  modes  (mechanisms)  .  For 
the  8-noded  Ahrad-Pawsey  element  these  are: 

1)  A  membrane  deformation  mode  with 


2)  A  bending  mode  with 

w  =  (1  -  ?2)  (1  -  n2>  -  j  U2  +  n2)  (75) 

Here  £  and  n  are  nondimens  iona 1  space  coordinates  chosen  such  that  they 
vanish  at  the  midpoint  of  the  element  and  equal  plus  or  minus  unity  on 
element  boundaries.  It  was  shown  by  Taylor  (see  Reference  6,  p.  189)  that 
the  bending  mode  is  not  "contagious”,  i.e.  it  occurs  only  for  a  single 
element.  With  more  than  one  element  the  global  stiffness  matrix  is  not 
singular  for  any  boundary  conditions  suppressing  the  menbrane  mode. 


Section  XV 


FLAT  ELEMENTS  FOR  CURVED  SHELL  ANALYSIS 


Curved  shell  elements  require  large  amounts  of  computer  time  for  formulation 
of  the  stiffness  matrix.  Also  the  use  of  curved  elements  introduces  the 
problem  with  self  straining.  Therefore,  curves  shells  have  most  frequently  been 
analyzed  by  use  of  flat  elements.  This  results  in  a  conformity  problem  that 
seems  to  be  particularly  important  in  nonlinear  or  stability  analysis. 

Figure  25  shows  two  flat  elements  representing  a  portion  of  cylindrical 
surface.  Rotational  compatibility  at  the  nodes  is  enforced  after  the  compon¬ 
ents  of  rotation  are  referred  to  the  same  set  of  coordinates,  implying  that 
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where  the  superscripts  (1)  and  (2)  refer  to  the  element  number. 

Mien  adjacent  flat  elements  meet  at  an  angle,  it  is  necessary  to  introduce 
the  normal  rotation  as  a  freedom  of  the  system.  A  disadvantage  with  this  is 
that  the  normal  rotation  does  not  appear  in  the  strain  energy  expression. 

As  a  consequence  the  resulting  linear  equation  system  becomes  increasingly 
ill-conditioned  as  the  angle  between  the  planes  of  adjacent  elements  becomes 
smaller.  Generally,  finite  element  computer  codes  define  a  small  limit 

a  and  if  a<a  ,  the  rotation  3,  is  ignored  and  as  an  approximation  the 

°  °  (1)  (2} 
conformity  constraint  becomes  3  =  3  . 


Figure  25.  Flat  Element  for  Curved  Shell  Analysis 


A  more  serious  problem  caused  by  the  use  of  flat  element  for  curved  surfaces 
is  related  to  interelement  displacement  nonconformity.  The  strain  energy 
expression  includes  derivatives  of  the  inplane  displacements  u,  v  of  first 
order  only  while  the  bending  strains  include  second  order  derivatives  of  the 
transverse  displacement.  Therefore,  w  is  usually  represented  by  polynomials 
of  higher  order  than  those  representing  u  and  v  .  Typically  w  is  cubic 
and  u  and  v  are  either  linear  or  quadratic  within  the  element. 


For  two  adjacent  flat  elements  at  an  angle  with  one  another,  complete  displace¬ 
ment  compatibility  requires  that  along  the  entire  boundary 


sin  (a/2)  =  0 

(77) 

(w(1)  -  w{2))  cos  (a/2)  +  (v(1)  +  v(2))  Sin  (a/2)  =  0 

Clearly,  these  conditions  cannot  be  satisfied  unless  v  and  w  along  the 
interface  are  represented  by  polynomials  of  the  same  order.  Failure  to  enforce 
this  compatibility  allows  the  individual  elements  to  buckle  under  axial  com¬ 
pression  as  plates  with  free  edges.  Displacement  conformity  can  be  restored 
by  introduction  of  constraints  on  w.  However,  this  results  in  a  very  stiff 
element  and  the  buckling  load  shows  rather  slow  convergence  from  above  with 
decreasing  grid  size. 


(v(1>  -  v<2’)  cos  (a/2)  -  (J»  *  „<2>) 


On  the  other  hand  displacement  conformity  can  be  restored  by  raising  the  order 
of  the  polynomials  representing  inplane  deformation.  With  a  cubic  representing 
w  it  is  necessary  that  third  order  polynomials  are  used  to  represent  u  in 
terms  of  the  y-coordinate  and  v  in  terms  of  the  x-ooordinate.  This  can  be 
achieved  through  inclusion  of  normal  rotation  ccrponents  as  supporting  freedoms 
for  the  inplane  displacement  field.  This  approach  was  used  in  the  development 
at  Lockheed  Missiles  &  Space  Company,  Inc  of  a  flat  element  for  thin  shell 
analysis.  The  out  of  plane  deformations  (bending  element)  were  defined  in  the 
same  way  as  in  the  QB12  element  discussed  above.  That  is,  rotational  nonconform¬ 
ity  is  permitted.  The  accuracy  is  of  second  order  for  rectangular  and  first 
ordc.  for  general  quadrilateral  elements. 


Inplane  displacements  were  defined  through  an  extension  of  the  bilinear  plane 
stress  element.  The  derivatives  u,x  and  v,  represent  normal  strains  and 
inclusion  of  these  as  degrees  of  freedom  would  lead  to  difficulties  in  the 
specification  of  boundary  conditions.  On  the  other  hand,  the  derivatives 
-v,x  and  u,y  represent  rotations  around  the  normal  of  short  line  segments  in 
the  x  and  y  directions,  respectively.  Unless  the  shear  strain  is  zero 
these  two  rotation  components  are  distinct  and  independent  of  one  another. 


The  shear  strain  is  equal  to  the  difference  between  these  rotation  ccrrponents, 
i.e. ,  u,  +  v,  .  If  the  shear  strain  is  allowed  to  be  disoontinous  at  the 

y  x 

node  there  can  be  more  than  two  distinct  normal  rotation  components.  Inclusion 
of  those  does  not  seem  to  be  practical  as  it  leads  to  great  complexities  in  the 
formulation.  After  introduction  of  the  normal  rotations  as  freedoms,  an  inplane 
displacement  field  is  obtained  in  which  displacements  normal  to  the  element 
boundaries  are  cubic  in  the  coordinate  along  the  boundary.  The  detailed 
derivation  of  the  stiffness  matrix  for  the  element  is  given  in  the  Appendix. 

The  use  of  the  normal  rotations  as  support  freedoms  for  the  inplane  displacement 
field  raises  the  order  of  the  polynomial  expressions  for  these  displacements  in 
one  of  the  coordinates  only.  As  the  constant  strain  elements  sure  known  to  have 
poor  convergence  properties,  the  orders  of  u  in  the  x-direction  and  v  in  the 
y-direction  were  raised  by  addition  of  tangential  displacement  freedoms  at 
midside  nodes.  The  element  is  referred  to  here  as  SH411.  Its  freedom  pattern 
is  illustrated  in  Figure  26. 

A  somewhat  simpler  version  of  this  element  SH410  excludes  the  midside  nodes  and 
uses  only  an  average  normal  rotation  component  at  the  comer  nodes.  This  means 
that  u  is  linear  in  the  x-direction  and  v  is  linear  in  the  y-direction,  and 
also  that  the  shear  strain  is  suppressed  at  the  nodes.  This  limits  the  useful¬ 
ness  of  the  element  to  special  cases.  The  SH411  element  has  nine  and  the 
SH410  has  five  integration  points.  However,  it  is  also  possible  to  use  reduced 
integration  with  four  points  in  each  of  these  elements. 

The  menbrane  part  of  these  elements  has  also  been  confcined  with  the  hybrid 
bending  element.  We  refer  to  these  elements  as  SH416  with  the  more  refined 
inplane  displacement  field  and  SH415  with  the  lower  approximation. 


The  problem  with  displacement  nonconformity  is  illustrated  here  by  use  of  a 
study  of  the  convergence  with  gridsize  of  the  buckling  loads  for  flat  and 
cylindrical  panels  subjected  to  axial  compression  or  shear.  Results  were 
obtained  with  the  elements  SH410  and  SH411  and  also  with  flat  elements  devel¬ 
oped  through  combination  of  the  bending  element  CF12,  discussed  above,  with 
quadrilateral  membrane  elements  obtained  through  combination  of  1ST  or  CST 
triangles.  We  refer  to  these  elements  as  SH420  (with  constant  strain 
triangles)  and  SH422.  FOr  flat  panels  the  buckling  mode  does  not  include 
inplane  displacements.  Therefore,  the  results  with  SH411  and  SH422  would  be 
practically  the  same  as  those  obtained  with  SH410  and  SH420  respectively,  and 
the  elements  with  higher  order  representation  of  inplane  displacements  are 
not  included  in  this  comparison. 

Results  of  an  analysis  of  buckling  of  a  flat  plate  under  axial  compression  are 
shown  in  Table  9 . 


TABLE  9 

BUCKLING  OF  FLAT  PLATE  IN  COMPRESSION 
Critical  Load  (lbs/in)  with 


Grid 

SH410 

SH420 

SH440 

3x3 

1391.1934 

(2.0) 

1454.2594 

(2.8) 

5x5 

1431.4386 

(4.2) 

1446.7181 

(7.6) 

1449.8566 

(8.1) 

7x7 

1439.5026 

(8.4) 

1446.2242 

(16.6) 

1442.9518 

(14.8) 

9x9 

1442.3723 

(20.1) 

1446.1372 

(29.8) 

1442.5526 

(25.5) 

11  x  11 

1443.7136 

(32.3) 

- 

1442.4792 

(40.4) 

13  x  13 

— 

- 

1442.4565 

(58.7) 

15  x  15 

1444.9247 

(69.9) 

- 

- 

19  x  19 

1446.5234 

(130.3) 

- 

- 

00 

1446.10 

- 

1442.42 

The  AHMAD- type  element  is  referred  to  here  as  SH440.  In  the  definition  of  a 
grid  the  mid  side  nodes  in  that  element  have  been  included,  i.e. ,  a  3  x  3  grid 
corresponds  to  only  one  SH440  element.  The  run  times  (CPH  on  the  CDC175 
NOS,  BE  system)  are  given  in  parenthesis  after  the  corresponding  buckling  load. 
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The  plate  dimensions  and  material  properties  are  the  sane  as  those  leading  to 
Equation  (74),  i.e.,  the  critical  load  is  1446.10  lbs/in.  With  the  SH440  the 
critical  load  is  somewhat  lower  since  transverse  shear  deformation  is  included. 
Extrapolation  from  the  results  in  Table  9  indicates  that  the  critical  load 
would  be  1442.42.  To  obtain  an  accuracy  in  the  budding  load  of  0.5  percent 
the  following  run  times  are  required 


with  SH410 

8.4 

sec 

SH420 

2.8 

sec 

SH440 

8.1 

sec 

SH415 

2.0 

sec 

If  a  higher-  accuracy  were  required  analysis  with  the  SH410  element  would  be 
considerably  more  expensive. 

Table  10  shows  results  from  budding  analysis  of  a  simply  supported  flat  plate 
in  shear.  The  plate  dimensions  are  5x5  inches  and  the  thickness  0.1  in.  and 
the  material  properties  are  E  =  10 ^  psi,  v  =  0.3.  The  reason  that  run  times 
are  smaller  in  this  case  in  comparison  to  those  with  the  same  grid  in  the  axial 
compression  case  is  that  the  convergence  in  the  eigenvalue  analysis  is  much 
faster.  The  results  with  SH415  indicate  fourth  order  accuracy  and  extrapolation 
from  these  values  we  determine  the  critical  loads  3371.3  lbs/ in.  without  and 
3339  lbs /in.  with  transverse  shear  effect  (SH440) .  The  figures  in  the  table  then 
shew  that  for  0.5  percent  accuracy  the  required  computer  time  is: 


with  SH410 

33 

sec 

SH420 

27 

sec 

SH440 

100 

sec 

SH415 

9 

sec 

If  the  accuracy  requirement  is  reduced  to  2  percent  the  corresponding  times  are: 


with  SH410 

11 

sec 

SH420 

15 

sec 

SH440 

40 

sec 

SH415 

6 

sec 

rinaa latent 
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TABI£  10 


BUCKLING  OF  FLAT  PLATE  IN  SHEAR 


Critical  load  (lbs/in)  with 


Grid 

SH410 

SH420 

SH440 

SH4L5 

5x5 

3182.59 

(3.6) 

3650.60 

(7.5) 

53063.0 

(7.3) 

3544.61 

(3.9) 

7x7 

3255.39 

(6.5) 

3423.80 

(15.5) 

— 

3406.72 

(7.1) 

9x9 

3301.26 

(10.9) 

3389.11 

(27.5) 

4093.58 

(19.6) 

3382.54 

(11.8) 

11  x  11 

3325.28 

(16.6) 

3379.00 

(44.1) 

— 

- 

13  x  13 

3338.92 

(24.1) 

— 

3403.48 

(40.2) 

- 

15  x  15 

3347.31 

(32.9) 

3373.30 

(80.9) 

3365.06 

(82.0) 

- 

17  x  17 

— 

3352.20 

(100.1) 

- 

19  x  19 

3356.61 

(52.7) 

3347.25 

(170.1) 

- 
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NG  OF  CYLINDRICAL  SHELL  SEGMENT  IN  COMPRESSION 


Results  for  a  cylindrical  shell  segment  under  axial  compression  are  shown  in 
Table  11.  The  segment  represents  a  complete  cylinder,  simply  supported  at  the 
edges,  with  a  length  of  11.0  in.  and  a  radius  of  36  in.  The  thickness  is 

7 

0.125  in..  Youngs  modulus  10  psi,  and  Poisson's  ratio  0.3.  It  is  assumed  that 
symmetry  conditions  prevail  at  midlength  and  that  the  cylinder  buckles  in 
15  circumferential  waves.  Consequently  a  5.5  in.  long  12  degrees  wide 
cylindrical  segment  is  considered  with  symmetry  enforced  on  three  sides.  The 
critical  load  according  to  an  analytical  solution  (Reference  24)  is  2500.69  lbs/in. 
Since  the  element  SH410  in  the  limit  supresses  inplane  shear  in  the  buckling 
mode  the  critical  load  converges  towards  a  value  that  is  somewhat  too  high. 

However,  the  results  are  surprisingly  close  to  the  analytical  solution  for 
all  grid-sizes,  the  error  in  no  case  being  above  0.3  percent. 

Both  for  SH420  and  SH410  it  appears  that  the  relatively  good  performance  is 
due  to  a  fortuitous  and  presumably  case  dependent  cancellation  of  errors  in 
different  directions.  Thus  it  makes  little  sense  to  include  those  in  a 
comparison  of  run  times  for  different  elements.  For  0.5  percent  accuracy  in 
the  buckling  load  the  following  run  times  are  required. 

with  SH411  13  sec 

SH422  not  achieved 
SH440  13  sec 
SH416  8  sec 

For  an  accuracy  of  2  percent  the  corresponding  results  are 

with  SH411  3  sec 

SH421  120  sec 
SH416  5  sec 

Since  the  elements  SH420  and  SH422  give  good  results  only  in  flat  plate 
analysis  and  in  sane  cases  may  perform  extremely  poorly,  they  are  not  included 
in  further  considerations.  However,  the  triangular  version  (the  Clough-Tocher 
triangle)  may  be  useful  in  shell  configurations  that  cannot  easily  be 
represented  by  quadrilateral  elements  only. 


In  order  to  include  effects  of  a  nonuniform  strain  distribution  in  the 
prebuckling  range,  a  case  was  also  considered  in  which  a  cylindrical  shell 
with  axisynmetric  imperfections  was  subjected  to  axial  compression.  The 
cylinders  radius  is  36  in.  and  its  thickness  0.125  in.  Material  data  are 
E  =  107  psi  and  \>  =  0^.  The  imperfection  is  defined  by 

Wq  =  0.03125  sin  (j  nx/L)  [cos  (tix/L)  -  1]  (78) 

Results  from  analysis  with  the  SH410  and  SH411  elements  are  compared  to 
results  with  finite  differences  (STAGSC)  in  Table  12 .  The  results  are  not 
quite  accurate  since  the  analysis  does  not  include  nonlinearities  in  the 
prebuckling  range.  However,  the  trend  is  quite  clear.  The  critical  load 
appears  to  be  about  1255  Ibs/in.  and  with  the  5x5  grid  and  SH411  elements 
a  2  percent  accuracy  is  obtained.  Similar  accuracy  with  the  other  two 
discretization  procedures  would  require  a  very  fine  grid  spacing,  in 
comparison  to  SH410  the  SH411  element  yields  the  same  accuracy  in  the  final 
results  with  about  one  fifth  of  the  run  time.  For  a  fixed  grid  size  the 
STAGSC  program  requires  much  less  computer  time  but  still,  the  run  time  is 
SH411  is  smaller  almost  by  an  order  of  magnitude. 

TABLE  12 

BUCKLING  OF  IMPERFECT  CYLINDRICAL  SHELL 

Critical  Load  (lbs/ in)  with 

Finite 

Grid  SH410  SH411  Differences 


3  x 

3 

— 

1651 

— 

4  x 

4 

1931 

1299 

— 

5  x 

5 

1578 

1280 

1865 

7  x 

7 

1405 

— 

1538 

8  x 

8 

1375 

- 

— 

13  x 

13 

1316 

— 

1346 
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Section  XVI 


INEXTENSIONAL  DEFORMATION 

Thin  shells  bend  easily  but  resist  stretching.  Therefore,  deformation  modes 
tend  to  be  inextensional.  For  a  straight  beam  element  the  neutral  surface 
strain  is 


£x  =  u'x  +  I  (u'x  +  w'x) 


(79) 


2  2 

Unless  the  rotation  of  the  element  is  very  large,  u,^  <<  w,^.  Inextensional 

deformation  then  implies  that  u,  is  approximately  equal  to  -  1/2  w,^.  With  a 

third  order  polynomial  representing  w  ,  w,  is  of  fourth  order.  In  that 

2  X 

case  the  relation  u,  =  -1/2  w,  can  hold  everywhere  within  the  element  only 
if  u  is  represented  by  a  polynomial  of  the  fifth  order.  This  problem  was 
first  recognized  by  Haftka  et  al.  in  Reference  28.  However,  accurate 
integration  of  the  nonlinear  terms  would  require  use  of  25  integration  points 
on  a  quadrilateral  plate  element.  This  leads  to  excessive  computer  time  and, 
since  reduced  integration  with  respect  to  the  nonlinear  terms  does  not 
introduce  mechanisms  in  the  system,  cubic  plate  bending  elements  are  usually 
based  on  a  nine  point  integration  scheme.  With  a  nine  point  integration 
scheme  it  is  possible  to  make  the  midsurface  strain  vanish  at  all  integration 
points  if  the  inplane  displacements  are  represented  by  cubic  shape  functions. 


The  simple  case  of  a  cantilever  beam  subjected  to  a  point  load,  as  shown  in 
Figure  27  can  be  used  for  demonstration  of  the  problem.  The  results  shown 
below  are  obtained  by  use  of  a  model  with  two  finite  elements.  The  lateral 
displacements  are  represented  by  a  cubic  and  a  three  point  Gaussian  scheme 
is  used  for  integration.  Linear,  quadratic  or  cubic  polynomials  are  used 
for  axial  displacements.  Within  the  limitations  of  the  theory  (a  Lagrangian 
formulation  and  the  curvature  defined  as  w,  )  the  case  with  cubic  u  gives 
exact  results.  The  case  of  linear  u  gives  poor  results  even  at  rather 
moderate  values  of  the  rotation  at  the  end  point.  A  two  to  three  degree 
rotation  appears  to  define  the  range  of  applicability  of  results  obtained 
with  such  an  element. 
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Figure  27.  Nonlinear  Bending  of  Beam 


The  results  for  a  beam  element  with  quadratic  representation  of  u  are 
close  to  the  accurate  solution.  However,  the  simple  case  considered  here  is 
not  quite  representative  and  analysis  of  more  ocrrplex  structures  show  that 
even  with  second  order  polynomials  for  u  and  v  ,  a  substantial  stiffening 
due  to  spurious  stretching  of  the  neutral  surface  can  take  place  even  at 
rather  moderate  rotations.  Ihe  element  SH411  discussed  in  the  preceding 
paragraph  can  be  modified  to  be  bicubic  in  the  inplane  displacement  components 
by  use  of  two  midside  nodes.  The  only  freedom  at  these  nodes  would  be  the 
tangential  displacements .  This  would  result  in  a  considerable  increase  in 
cemputer  time  and  it  appears  that  the  use  of  reduced  integration  may  be  a 
more  attractive  solution. 

Reduced  integration,  with  four  Gaussian  points  (i.e. ,  reduction  beyond 
exact  integration  of  the  second  order  terms  in  the  energy)  was  introduced  as 
an  option  in  the  elements  SH410  and  SH411.  Figures  28  and  29  show  some 
results  on  nonlinear  bending  of  a  thin  and  relatively  deep  arch.  The  thick¬ 
ness  is  0.125  in.  and  the  width  1.0  in.  A  solution  to  this  problem, 
presented  in  Reference  29,  is  indicated  by  the  dotted  curve.  The  results 
with  SH411  elements  are  also  shown  in  Figure  28.  Only  half  the  arch  is  modeled 
and  the  analysis  was  performed  with  either  four  or  eight  elements  on  this 
model.  Clearly  the  rate  of  convergence  with  grid  size  is  much  better  if  the 
reduced  integration  is  used,  so  that  spurious  stretching  of  the  middle 
surface  is  avoided. 

Results  from  analyses  with  SH410  elements  are  shewn  in  Figure  29.  The 
difference  between  use  of  four  or  five  integration  points  is  surprisingly 
large  but  even  with  four  points  spurious  neutral  surface  stretching  cannot 
be  avoided  and  the  results  seem  to  have  little  resemblance  to  actual 
structural  behavior. 

Although  significantly  different  the  results  obtained  here  as  well  as  those 
in  Reference  29  appear  to  represent  converged  solutions  (with  grid  size) . 

Most  likely  the  discrepancy  is  due  to  differences  in  the  basic  theory.  For 
example  the  definition  of  the  strain  in  Reference  29  does  not  include 


Figure  29.  Load  Displacement  Curves  for  Arch,  SH410  Elements 
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the  term  1/2 (v,^)  .  The  formulation  for  SH411  w?is  tenporarily  modified  through 
exclusion  of  this  term.  With  the  modified  elemea.ts  the  computed  limit  point 
is  only  two  percent  belcw  the  one  reported  in  Reference  29.  The  deformed 
shape  at  the  limit  point  as  computed  with  SH411  (4  elements,  4  integration 
points)  is  shown  in  Figure  30. 

Further  exploration  of  the  use  of  reduced  integration  to  avoid  spurious  stretch¬ 
ing  of  the  middle  surface  indicated  even  greater  advantage  for  more  complex 
cases.  For  example  a  good  estimate  of  the  collapse  load  for  long  cylinders 
in  bending  (Brazier  effect)  cam  be  obtained  with  a  relatively  coarse  grid 
(6x6  elements  for  the  exanple  used)  with  SH411  and  with  four  integration 
points.  With  nine  integration  points  and  the  same  grid  the  collapse  load 
is  overestimated  by  a  factor  of  two. 

An  attempt  was  also  made  to  solve  the  arch  problem  (Figure  30)  by  vise  of  the 
Ahmad  type  element,  SH440.  However,  even  with  a  very  fine  grid,  21  points 
over  the  half  arch  the  load  displacement  curve  does  not  contain  a  limit  point. 

It  seems  that  this  element  has  a  tendency  to  "lock"  with  increasing  rotation 
of  structural  elements.  At  the  limit  point  the  maximum  rotation  in  the 
deformation  pattern  is  about  35  degrees.  Possibly  the  locking  of  the  element 
is  due  to  the  fact  that  the  transformation  of  a  vector  of  rotation  components 
is  valid  only  in  the  small  rotation  range.  However,  displacement  configura¬ 
tions  with  rotations  as  large  as  35  degrees  are  not  within  the  range  of  the 
moderate  rotation  theory.  In  that  case  all  elements  using  rotations  as 
freedom  would  suffer  from  this  problem  in  large  displacement  analysis  based 
on  a  Lagrangian  formulation.  It  appears  that  the  elements  SH410  and  SH420 
elements  are  less  sensitive  to  the  locking  problem.  More  light  is  thrown 
on  this  problem  in  the  following  section. 
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Section  XVII 


SHELL  Cm  APSE  ANALYSIS 


For  a  final  evaluation  of  the  shell  elements,  represented  by  SH411,  and  the 
degenerate  brick  element  (SH440)  a  couple  of  cases  with  strongly  nonlinear 
behavior  were  selected.  The  first  one  of  these  is  the  so  called  "pear-shaped” 
cylinder,  first  presented  as  an  example  case  in  Reference  30.  load  displace¬ 
ment  curves  (Figure  31)  cure  shown  from  analysis  with  the  SH411  elanent  for  three 
different  grid  configurations.  The  conputed  failure  loads  are 


with  3  x  27  grid  3586  lbs 
5  x  37  grid  2731  lbs 
7  x  47  grid  2586  lbs 


For  moderate  values  of  the  load  there  is  no  significant  difference  between 
the  results  for  different  grid  sizes.  The  convergence  difficulties  at  higher 
load  levels  may  be  due  to  the  fact  that  the  displaconent  pattern  becomes 
more  complex  with  higher  values  of  the  load.  It  could  also  be  related  to  the 
locking  phenomenon  discussed  in  the  previous  section. 

load  displacement  curves  for  the  pear-shaped  cylinders  obtained  by  use  of  the 
SH440  element  are  shewn  in  Figure  32.  The  failure  loads  cure 


with  5  x  37  grid  2657  lbs 
7  x  47  grid  2530  lbs 

Since  the  elements  appear  to  have  a  tendency  to  lock  with  increasing  values 
of  the  rotation  components,  the  problem  of  convergence  with  grid  size  becomes 
complex.  It  is  not  possible  to  relate  the  error  to  the  grid  size  alone  and 
extrapolate  with  confidence.  However,  it  appears  that  in  both  cases  the 
collapse  load  is  somewhere  around  2300  or  2400  lbs,  that  is  close  to  the 
2340  lbs  obtained  with  a  STAGS  A  analysis  based  on  a  9  x  45  grid  reported 
in  Reference  30.  With  the  same  grid  size  the  finite  difference  program 
gives  results  that  are  at  least  as  close  as  those  obtained  by  use  of 


AO  (lb) 


Figure  31.  Load  Displacement  Curves  for  "Pear-Shaped"  Cylinder  with 
SH411 
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AXIAL  LOAD  (lb) 


Figure  32.  Load  Displacement  Curves  for  "Pear-Shaped"  Cylinder  with 
SH440 
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any  of  the  two  finite  elanents,  at  a  cost  that  is  an  order  of  magnitude  less. 

The  second  case  considered  is  the  cylinder  with  a  rectangular  cutout  for  which 
analytical  as  well  as  experimental  results  are  presented  in  Reference  31.  The 
dimensions  of  the  cylinder  and  the  results  from  analysis  with  the  SH411  element 
are  shown  in  Figure  33. 

The  computed  failure  loads  are 

with  11  x  13  grid  3698  lbs 
17  x  19  grid  3060  lbs 
23  x  27  grid  2922  lbs 
29  x  37  grid  2750  lbs 

With  SH440  elements  the  following  results  were  obtained 

17  x  19  grid  3733  lbs 
23  x  29  grid  3098  lbs 

Finite  difference  solutions  and  experimental  results  (Reference  31)  are  in 
good  agreement.  A  finite  difference  analysis  with  a  21  x  23  grid  gives  a 
collapse  load  of  £250  lbs.  in  good  agreement  with  experiment.  It  is  not  econ¬ 
omically  feasible  to  obtain  a  converged  solution  in  this  case  with  any  of  the 
finite  elements .  The  true  value  of  the  critical  load  may  be  somewhat  higher 
than  that  obtained  in  the  experiments,  but  still  it  appears  that  the  finite 
difference  approach  is  vastly  more  efficient  for  this  type  of  analysis .  The 
maxinun  rotations  at  failure  for  the  two  cases  discussed  in  this  paragraph  are 
rather  moderate ,  about  7  degrees  in  both  cases . 
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LATERAL  DISPLACEMENT  AT  CUTOUT  EDGE  (In.) 


Figure  33.  Lead  Displacement  Curves  for  Cylinder  with  Cutout  with 
SH411 


96 


Section  XVIII 


CONCLUSIONS  AND  REQCtMENDKTIONS 


More  than  anything  else  the  results  presented  in  the  preceding  sections  tend 
to  stress  the  importance  of  the  choice  of  discretizatijn  procedure  in  analysis 
of  shells,  particularly  when  it  comes  to  buckling  and  collapse  analysis.  It 
is  shewn  for  example  that  the  price  for  a  buckling  analysis  with  a  two  percent 
accuracy  of  a  cylindrical  panel  may  vary  by  as  much  as  a  factor  of  40 
depending  on  the  choice  of  element.  While  analysis  with  SH410  requires  about 
half  the  run  time  in  comparison  with  SH411  for  buckling  of  a  panel  with  uniform 
prestress,  it  requires  about  five  times  as  much  if  an  axially  syrrmetric 
inperfection  is  introduced  to  generate  a  nonuniform  prebuckling  stress  distri¬ 
bution.  Most  striking  may  be  is  that  while  the  finite  element  configurations 
investigated  here  result  in  very  efficient  analysis  of  bifurcation  buckling 
loads  they  appear  to  be  hopelessly  inferior  to  the  finite  difference  method 
(elements  with  extended  support)  in  nonlinear  analysis  of  shells  that  undergo 
relatively  large  rotations  (5  to  10  degrees)  before  collapse  occurs.  It  seems 
reasonable  to  assume  that  this  problem  is  related  to  the  use  of  shell  equations, 
based  on  the  moderate  rotation  assumption  together  with  a  Lagrangian  formulation. 
The  finite  difference  procedure  does  not  include  rotations  as  freedoms.  This 
may  explain  why  it  is  less  sensitive  to  these  assumptions . 

For  bifurcation  buckling  analysis  of  thin  shells  the  elements  SH415  and  SH416 
are  efficient  and  reliable  (SH416  only  if  the  prebuckling  membrane  strain  is 
reasonably  uniform) .  While  sometimes  the  elements  SH410  and  SH411  may  lead  to 
semevhat  more  economic  analysis,  due  to  the  balance  between  positive  and  neg¬ 
ative  errors,  they  are  for  the  same  reason  less  reliable.  In  particular  the 
tendency  to  converge  from  below  in  many  cases  is  undesirable  as  it  can  lead  to 
spurious  local  buckling  behavior  unless  a  very  fine  grid  is  used.  Also  SH410 
and  SH411  are  not  suitable  if  the  element  planform  deviates  much  from  a  rectangle. 
On  the  other  hand  the  hybrid  forms  SH415  and  SH416  can  not  be  used  without  some 
penalty  on  computer  time  or  accuracy  for  shells  exhibiting  coupling  between 
membrane  and  bending  behavior.  For  such  cases  the  elements  SH410  and  SH411  are 
recomnended  if  the  elements  are  close  to  rectangular. 

It  must  be  stressed  that  the  investigation  reported  here,  although  extensive  ,  is 
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far  from  exhaustive.  For  exanple  Iran's  Semi-loof  element  (Reference  32)  and 
the  hierarchy  of  modified  degenerate  brick  elements  recently  developed  by 
Hughes  and  his  coworkers  (Reference  33)  should  be  considered.  No  curved 
elements  based  cn  shell  theory  were  included  in  the  study.  Except  for  flat 
plates  then,  the  analysis  based  on  the  shell  elements  included  is  reduced  to 
a  second  order  accuracy,  that  is  the  same  as  in  the  finite  difference  analysis. 

For  thicker  shells  the  analysis  must  be  based  on  a  second  order  theory.  For 
this  case  the  Ahmad  type  element,  SH440,  is  the  only  element  included  in  the 
investigation  that  will  lead  to  accurate  results  (for  bifurcation  budding  as 
well  as  nonlinear  collapse  analysis).  Generally,  the  Ahmad  element  works  well 
also  for  thin  shells,  and  the  fact  that  it  is  the  only  choice  in  seme  shell 
codes  has  same  justification.  It  is  feasible  that  the  relatively  slow  conver¬ 
gence  in  plate  shear  buckling  will  be  cured  by  inclusion  of  the  ninth  node  (at 
the  midpoint).  However,  with  reduced  integration  this  may  lead  to  problems 
with  mechanisms.  Also,  it  may  be  noted  that  the  Ahmad  type  elements  will  be 
penalized  in  terms  of  computer  time  for  shell  walls  with  membrane-bending 
coupling. 

The  most  disturbing  aspect  of  this  investigation  is  undoubtedly  the  poor  per¬ 
formance  of  the  finite  elements  included  in  the  collapse  analysis  of  the  pear- 
shaped  cylinder  and  the  cylinder  with  a  cutout.  Since  the  finite  elanent  for¬ 
mulation  with  the  same  grid  size  gives  a  more  accurate  description  of  the  de¬ 
formation  pattern  it  should  not  require  a  finer  grid  in  the  collapse  analysis. 
There  are  good  grounds  therefore  to  assure  that  the  apparent  stiffening  of  the 
structure  with  increasing  rotation  can  be  cured  by  a  nonlinear  analysis  procedure 
that  includes  "updating  of  the  geometry".  The  so-called  "corotational  approach" 
advocated  by  Horrigmoe  and  others  (Reference  34,  for  exanple)  appears  to  be  most 
promising  since  it  allows  for  such  updating  without  modification  (with  increasing 
deformation)  of  element  shape  functions. 


It  is  clear  that  much  work  remains  to  be  done  in  the  area  of  discretization 
procedures  for  stability  and  nonlinear  collapse  analysis  of  thin  shell  structures 
It  is  recommended  here  that  flat  and  curved  shell  elements  as  well  as  degenerate 
brick  elements  be  introduced  in  STAGSC-1  for  analysis  of  an  extended  set  of 
benchmark  cases .  With  the  corotational  approach  the  problem  with  self-straining 
of  curved  elements  should  be  essentially  eliminated  while  a  third  or  fourth 
order  accuracy  is  maintained.  If  such  elements  do  not  prove  to  be  substantially 


98 


t 

I 


r 


1 


more  economic  than  the  flat  elements  for  curved  shell  analysis  the  possibilities 
of  reintroducing  finite  difference  formulations  must  be  seriously  considered. 
However,  this  alternative  seems  less  attractive  because  much  work  vrould  be 
needed  in  order  to  derive  formulations  that  are  comparable  to  finite  elements 
in  modeling  capability.  Particularly,  programming  procedures  involving  geometry 
updating  for  problems  with  large  rotations  will  be  considerably  more  complex 
if  the  elements  are  not  self-contained. 
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Appendix 

STIFFNESS  MATRIX  FOR  THE  OUAF  ELEMENT 


The  degrees  of  freedom  of  the  system  are  the  displacement  vectors  u8, 
s  s 

v  ,  and  w  on  the  shell  (or  plate)  and  the  vector  of  rotation  components 
£  £  _ 

Y1  '  y2  *  ®  1  an<*  ®  2*  aS  s^own  ^*8ur®  Al.  These  vectors  represent  the 
components  of  displacements  and  rotations  in  a  local  coordinate  x',  y',  z' 
system  where  the  x'  and  y'  axes  are  tangents  to  the  coordinate  lines,  and 
z'  Is  directed  along  the  outward  normal.  The  component  of  lnplane  rota- 

8  S 

tlon,  or  Y2  »  through  the  node  point.  If  a  set  of  surface  coordinates 

£ 

x,  y  are  used  for  the  grid  generation,  the  notation  y^  refers  to  the 
rotation  of  a  coordinate  line  corresponding  to  constant  y  and  y^s  refers 
to  the  rotation  of  a  line  corresponding  to  constant  x. 


Fig.  Al  The 
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Next,  a  flat  aleaent  la  considered  with  the  corners  at  four  node  points 
as  shown  in  Figure  A2.  A  total  of  eight  nodes  are  defined  on  the  element 
as  shown  in  the  figure.  A  Cartesian  system  x,y,z  is  Introduced  with  x,y 
.in  the  directions  as  shown,  and  z  completing  a  right-handed  system. 


[7]  -  Side  NuM&EfL  (7)  -  Node  Nm&ce. 

Fig.  A2  The  Flat  Element 


The  3x3  matrix  transforms  a  vector  in  the  system  x",  y',  z'  at  corner 

i  into  a  vector  in  the  coordinate  system  x,y,z  for  the  element.  The 
displacement  components  at  each  of  the  corners  are  obtained  from 


tt) 

f"8] 

- 

.  [T3,i) 

v8  1 
1  ( 

) 

s 

V  wj 

(Al) 
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corresponds  to  the  rotetlon  of  side  number  i  at  corner 


The  vector 
number  j .  Thus 


Pi 

h 

e 3 


(3.  j) 


} 


^2 


0) 


j  =  2.3 


(A2) 


j  =  3,4 


The  next  atep  la  to  determine  the  derivatives  of  the  displacement  compon¬ 
ents  with  respect  to  the  coordinates  €  and  n  which  are  defined  so  that 
K  •  -1  on  side  4  and  4-1  for  aide  2  and  n  *  -1  on  side  1  and  +1  for  side  3. 


A- 3 


•  >-  ■ 


That  is 


x  *  <0^>  {x^} 

y  =  <0*>  fy) 

where 

<0(1)>  =  <i/4(i-5)(i-<n),  i/4(i+§(i--n), 

1/4(1+5)(1+T1)»  1/4(1- §)(1+T))> 


(A3) 


where  x^,  are  the  x  and  y  coordinates  of  the  four  corners.  First  the 
derivatives  are  determined  with  respect  to  a  parameter  representing  the 
distance  along  the  side  of  the  element,  see  Figure  A3. 


-  y,v 


Fig.  A3  Inplane  Rotation  of  x-Coordinate  Line 


For  example,  for  line  1  at  corner  1 


g!l)  =  p/1, 2)  sin  «,  +  B 


+  »  O.D 

1  +  P2  COS  Q  j 


(A4) 
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With  the  notations 


c  = 

cos  (a^)  and 

s  = 

•i 

ain  (a^)  at  corner  1: 

II 

A3|(t> 

2(x2  -x^/c^i 

and  consequently 

3w<l>  , 

Z(x2-xi) 

s  *«*•»♦  e  P,(1* l) 
a^  1  a^  2 

o F 

cal 

The  derivatives  with  respect  to  n  are  determined  in  the  same  way  (See 
Figure  A4). 


Fig.  A4  Inplene  Rotation  of  y-Coordlnate  Line 


Corresponding  derivative!  at  the  remaining  three  corners  ere  then  deter-* 
■lned  In  the  ease  way. 

The  field  of  lateral  displacement  Is  determined  from  the  lateral  displace¬ 
ments  and  their  derivatives  with  respect  to  the  Cartesian  coordinates  x 
and  y  at  the  corner  nodes.  The  chain  rule  of  differentiation  yields 


The  transformation  matrix  J*  is  obtained  as  the  transpose  of  the  inverse  of 


[J] 
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and 


[J*]  = 


ay  ax 

Vn  ’  a  7i 


ay. 

as 


ax 

aTJ 
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ay  ax 

an  ag 


dy. 

as 
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The  derivatives  of  x  and  y  with  respect  to  C  and  n  are  readily  obtained 
from  Eqs.  (A3).  The  freedoms  of  the  flat  element  (in  addition  to  v  and 
its  derivatives)  are  u  and  v  displacements  at  the  corner  nodes,  the  cor¬ 
rection  to  the  inplane  tangential  displacements  at  nodes  5,  6,  7,  and  8, 
and  the  corner  node  freedoms  corresponding  to  the  boundary  line  rotations 
about  the  normal  to  the  plate  (see  Figure  Al) . 


The  displacement  field  that  matches  these  freedoms  will  now  be  defined. 

A  bilinear  field  is  used  to  match  the  corner  displacement.  The  biquad¬ 
ratic  displacement  field  has  zero  displacements  at  corner  nodes.  It 
represents  an  addition  to  the  tangential  displacements  at  the  midsize 
nodes  as  determined  from  the  bilinear  field.  The  parameters  yu  and  yv 
are  the  corrections  to  the  displacement  field  that  are  necessary  in  order 
that  the  inplane  rotations  will  match  the  rotations  of  boundary  lines  at 
the  corners  (as  determined  byg^,  see  Eqs.  A2). 


Hence,  we  have 


u  = 


v  * 


<0(1)>  fu.}  + 
<0(1)>  fv.)  + 


<0(2)>  (t.uJ  +  <0(3)>{y“}  +  <0(3)> 

1  l  *  U  1  ’l  *  vu 

<0(2)>  ft  V)  +  <  0(3)>  fvVl  +  <  0(3)> 
i  }  v  'i  '  uv 


i  =  1,  2,  3,  4  (A15) 

Notice  that  each  of  the  four  costponents  of  the  biquadratic  field  corres¬ 
ponds  to  displacements  parallel  to  one  of  the  aides.  Likewise,  the  com¬ 
ponents  of  the  bicubic  field  include  only  displacements  in  a  direction 
normal  to  one  of  the  sides.  Consequently 
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(A16) 


ftU}  =  <t.  c  ,  t,  a  ,  t.  c  ,  tA  b  > 

1  1  J  1  2  »2  3  o  4  of^ 

ft.Vl  =  <  -t,  ■  ,  t,  C  ,  -t,  B  ,  t.  C  > 

li  1  1  ofj  2  or2  3  ofj  4 

»  The  form  of  the  components  in  the  bilinear  field  are  given  by  <0^>  in 
Eqs.  A3.  That  is, 

<0(2)>=  <1/2(1-52)(1-T]).  1/2<1+SK1-712),  1/2(1- SZ)<1+T|).  1/ 2 (1  - 9 (1  - T]2)  > 

(A  17) 

and 

<0(^>  =  <1/16  ca4  (2-3C+53)  (1  -Tl-Tl2+n3),  1/16  ca2(2+3!-F3)(l-7]-T|2+T)3), 

1/16  ca2  (2+3?-53)(-l-T)+T)2+Tl3).l/l6Ca4(2-3E+?3)(-l-W+'n3)> 

<0Jtv>=  <-1/16  *cr4U-3?+?3)  (l-TlV+Tl3),  -1/16  So2{2+3§-?3)(1-T)-T12+T13), 

-1/16  sq2  (2+ 3?-£3)  (-1  -T)+7)2+T)3),  -1/16  sq4(2-3§+§3)  (- 1 -T)+t2+713)  > 

<0^3,>=  <1/16  col  (l-?+?2+S3)  (2-3D+7)3).  1/16  C(yl(-l-5+|2+53)(2-3T)+713), 

1/16  c^3(-l-|+?2+e3)(2+3Tl-Tl3)>  1/16  Ctt3(l-?-?2+^3)(2+3n-T13)> 

<0^  >  -  <l/16«wl  (l-?+i;2+?3)(2+3Tl-Tl3).  1/16  stt3(l-P-e2+F3)(2+3TrTl2)> 

(A  18) 

Differentiating  Eq.  (A15)  with  respect  to  the  normalized  coordinate  yields 
in  matrix  notation 


where 


T  =  1/2 
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,  T  =1/2 
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-1 
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-1 
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If  Che  operations  Indicated  by  the  second  tern  on  the  right  side  of  Eq. 
(A19)  are  carried  out,  this  tern  is  found  to  be  identically  equal  to 
zero  at  all  corner  nodes.  Consequently,  the  derivatives  with  respect  to 
€  and  n  of  the  displacement  field  corresponding  to  the  midnode  displace¬ 
ments  will  vanish  and  need  not  be  given  any  further  consideration. 


Next  the  boundary  line  normal  rotations  corresponding  to  the  linear 
polynomial.  For  this  purpose  the  Inplane  displacements  normal  to  the 
boundary  are  differentiated  with  respect  to  the  distance  along  the  ele¬ 
ment  side. 


From  the  figure,  it  la  observed  that 
r  =  S(x2-xt)/2  . 

8  =  T)(y4-yi)/2  .  ^ 


2  jj_ 

(x2-xL)  95 


2  _a_ 

(y4  -  yj) 


on  boundary  1 

(A21) 

on  boundary  4 


Also, 


v 

n 


v  c 


+ 


+ 


u  s 


*1 
u  c 

a4 


on  boundary  1 


on  boundary  4 
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From  comparison  of  the  left-hand  side  of  Eq.  (A25)  to  similar  terms  In 
Eq.  (A19)  the  expression  for  the  correction  terms  yU  and  yv  are  deter¬ 


mined,  yielding 


where  j  =  1,  if  i  =  1  or  2,  j  =  3  if  i  =  3or  4,  and  k  =  4  if  i  =  1  or  4,  k  =  2  if 
i  =  2  or  3. 


Tuu*  Tw’  V  TPv  ar*  defined  in  Etls-  <A2°).  (A2&).  <A27);  the  j>  are 
defined  in  Eqs.  (A2).  Notice  also  that  u,g  =  -63^  °r  **  ^  and  v,f  =  83^  °r 


(A31) 


By  substituting  Eq.  (A30)  into  Eq.  (A13),  the  correction  terms  y11  and  yV 


are  eliminated  and  the  expression  for  u  and  v  Is  obtained  in  terms  of 
the  selected  flat  plate  freedoms. 


A-12 


„  .  <«<l)>f  u,l  ♦  <*(2)>  ft,"} 


v  = 


♦  <tW>  <Tuuf“l1tTvurv.)*TBu!u.Sil) 

+  <*w>  + 

(Vjj  +  >  [t/) 

+  <»l3'>  lT»v'Uil+Twtvi’+Tevfv'ri’> 

+  <>v’>  cr„f,)  +  T„tv,MTerCv.r.l> 


The  derivatives  of  u  and  v  with  respect  to  x  and  y  are  needed  at  the 
Integration  points  inside  the  element.  The  chain  rule  of  differentia¬ 
tion  yields 


v,. 


J* 


J* 


where  J*  is  defined  by  Eq.  (A14). 
Hence, 


J  =(1/K) 


l'Tl 
'*  Tl 


0>>.H{yi> 

■  «  *">.  5  <yi> 

A 

V 

• 

r~% 
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(A  32) 


(A33) 


(A34) 


w-  Y4i 


where 


k  *  <<#U)>tn  (yl))  << 0<l)>  5 u.] > 

-(<  0(1)>,  l  ^i>)  (<  0<1>>,  Tt  <A35) 

The  values  of  the  shape  functions  0^,  and  0^  at  the  eight  nodes 

(four  corners  and  four  mid  nodes)  are  shown  in  Table  Al. 
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Table  A1 

SHAPE  FUNCTIONS  AND  DERIVATIVES  EVALUATED 
AT  NODES  AND  MID  NODES 

a)  Bilinear  Functions 

=  <-i/4(i--n),  i/4(i-m,  i/^d+Ti).  -i/4o+n)> 

0  5 

*  <-l/4(l-5).  -l/4(l*«.  1/4(1H).  l/4(l-5)  > 

•<v  >,11 


Node 

5 

1) 

•w 

»!"t 

,1\ 

1 

-1 

-l 

1,0, 0,0 

-1/2, 1/2,  0,  0 

-1/2,  0,0, 1/2 

2 

1 

-1 

0.1.  0,0 

-1/2, 1/2,  0,0 

0, -1/2, 1/2. 0 

3 

1 

l 

0.0,1.  o 

0,0, 1/2,  -1/2 

0,  -1/2, 1/2,0 

4 

-1 

1 

0,0, 0,1 

0, 0,1/2, -1/2 

-1/2,  0,0. 1/2 

5 

0 

-1 

1/2. 1/2,0, 0 

-1/2, 1/2,  0,0 

-1/4,  -1/4, 1/4,  1/4 

6 

1 

0 

0,1/2, 1/2,0 

-1/4, 1/4, 1/4,  -1/4 

0,  -1/2, 1/2,0 

7 

0 

1 

0,0, 1/2, 1/2 

0,  0, 1/2-1/ 2 

-1/4,  -1/4, 1/4, 1/4 

8 

-1 

0 

1/2.0, 0,1/2 

-1/4, 1/4, 1/4,  -1/4 

-1/2,  0,0. 1/2 
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b)  Biquadratic  Function 


<0(2)>  =  <-§(l-TU.  1/2(1-T12).  -g(l+tl).  -1/2(1-T12)  > 

•  ? 

<0<2)>  =  <-l/2(l-?2),  -(Mm  1/2(1'52).  -(1-5)  > 

»  M 


Node 

5 

1) 

*<2> 

0(2) 

\5 

0<2> 

.  T| 

1 

-1 

-1 

0,0,0,  0 

2,  0,0,0 

0,0,0,  2 

2 

1 

-1 

0,0, 0,0 

-2,  0,0,0 

O 

o 

w 

o 

i  3 

1 

1 

1 

0,0,  0,0 

0,0.  -2,0 

0, -2.0,0 

! 

4 

-1 

1 

0,  0,0,0 

0,0.  2,0 

0,0,  0,  -2 

'  5 

0 

-1 

1,0,  0,0 

0,0,  0,0 

-1/2,1, 1/2,1 

i  <■ 

l 

0 

0,1, 0,0 

-1,1/2,  -1.  -1/2 

o 

o 

o 

o 

1 

1  7 

i 

0 

1 

0,0, 1,0 

0,0,  0,0 

-1/2,  -1,1/2,  -1 

8 

-1 

0 

0,0, 0,1 

1, 1/2,1,  -1/2 

o 

o 

o 

o 

c)  Bicubic  Functions 


<  *<3)>  , 
u  ,r\ 


<1/16  c  (-3+3S2)(l-H-T|2+Tl3),  1/16  c  (3-3g2)(l-*n-Tf+lf). 

®4  °Z 

1/16  c  {S-SS^t-l-D+^+T?),  1/16  c  (-3+3C2H-l-Tl+‘n2+Tl3)> 
*2  ®4 

<1/16  c  (2-3§+S3)(-l-2T)+3T)2),  1/16  c  (2+ 3 g - g3) { -1-2 T)+3 T)2), 
®4  aZ 

1/16  c  (2+3g-g3)(-l+2T\+3T)2},  1/16  c  (2-3 g+ g3) (-1+2T1+3 T\2) 
aZ  ®4 

<1/16  c  <-l-2g+3^) (2 -3T.+T?),  1/16  c  (-1+2 g+3 g2) (2 - 3 7)+ Tl3) 

®1  *1 

1/16  c  (-l+25+3S2)(2+3H-T\3),  1/16  c  (-l-2?+3§2)(2+3Tl-‘n3)  > 
®3  “3 


<0-3)>.H  *  <1/16  C~  (1-5-52+5JH-3+3^2).  1/16  c„  (-l-g+g2+g3)(-3+3-n2). 


a-j  ' 

1/16  C/y  (-l-g+g2+g3)(3-3Tl2),  1/16  c  (l-5-g2+g3)(3-3Tl2)  > 
®3  *3 


<0 


(3) 

uv 


<0<3> 


> 


*  <-i/i6  §  (-3+3g2)(i-*n-*n2+-n3).  -1/16 .  (3-3?2Hi-Ti-if+T?), 

®4  ®2 

-1/16  s  (3-352)(-I-H+Tl2+n3),  -1/16  •  (-3+§2K-l-Tl+Tt2+713)> 

“2  ®4 

*  <-1/16  •  <2-3g+g3)(-l-2Tl+3Tl2),  -1/16  (2+3g-g3) (-1-2T1+3 Tl2). 

®4  ®2 

-1/16  s  (2+3g- g3)(-l+2T>+3*n2).  -1/16  s  (2 -3g+ g3) ( -1+2T1+3T12)> 

®2  aZ 

=  <1/16  s  (-l-2g+3g2)(2-3T)+T|3),  1/16  s  (-l+2g+3g2)U~3T]+T)3). 

®1  ®1 

1/16  s  3  ( -l+2g+3  g2) (2+371-TJ3) ,  1/16  s  < -1 -2g+3  g2)  (2 +3  Tl-*n3)  > 

o  a  j 

*  <1/16  •  (l-g-52+g3)(-3+3*n2).  1/16  •  <-l-g+g2+g3)(-3+3712). 

*1  ®1 

1/16  (-l-5+52+^)(3-3H2).  1/16  s  (l-g-g2+g3)(  3-3Tl2)> 

*3  *3 
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The  changes  of  curvature  in  the  element  are  determined  from  the  lateral  dis¬ 
placements  and  the  rotation  components  w,  and  v,  as  given  by  Eq.  (All). 

x  y 

The  cubic  interpolating  functions  defined  for  the  membrane  element  could  be 

•used  If  the  analysis  were  restricted  to  rectangular  elements.  However,  for 

the  bending  part  the  Isoparametric  mapping  to  the  C,n  coordinates  produces 

strain  energy  under  rigid  body  displacement  if  the  element  is  nonrectangular. 

Therefore  the  Taylor  series  approach  was  used.  The  two-dimensional  Taylor 

series  includes  all  terms  through  the  third-order  and  two  fourth-order 
3  3 

terms  w?^x  y  and  v^xy  .  A  penalty  function  is  applied  to  the  freedoms 
wi;j  and  w^.  This  constraint  is  found  to  reduce  significantly  the  effect 
of  rotational  nonconformity. 


